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CHAPTER  1 
INTRODUCTION 

Advanced  f fber-reinforced  composite  materials  such  as 
boron/epoxy  and  graphi te/epoxy  have  been  successfully 
employed  as  structural  materials  In  al rcraf ts , ,  ml ss 1 1 es  and 
space  vehicles  In  recent  years,  and  the  performance  of  these 
composites  has  shown  their  superiority  over  metals  in 
applications  requiring  high  strength,  high  stiffness  as  wel 1 
as  iow  weight.  The  advantages  of  these  compos  I tes ,  however, 
are  overshadowed  by  their  relatively  poor  resistance  to  the 
impact  loadings,  which  has  prevented  the  application  of 
these  materials  to  turbine  fan  bladings.  Many  other  reports 
deaung  with  the  responses  of  advanced  composites  to  various 
typos  of  impact  have  further  increased  the  need  for  a  better 
understanding  of  the  problem  so  that  the  survivability  of 
these  composites  can  be  improved. 

It  is  obvious  that  impact  produces  damage  and 
consequently  reduces  the  strength  of  composite  materials. 
The  damage  modes  usually  include  local  permanent 
deformations,  breakage  of  fibers,  del  ami nat ions ,  etc.. 
While  the  cause  of  these  damages  are  still  unknown  and  may 
not  be  simple  in  nature,  in  general,  the  impact  of  a  soft 
object  could  give  a  longer  contact  duration,  and  the  dynamic 


response  of  the  whole  structure  is  of  importance.  The  hard 
object  impact  usually  gives  a  short  contact  time  and  results 
In  the  initial  transmisson  of  impact  energy  Into  a  local 
region  of  the  structure.  This  initial  energy  will  propagate 
into  the  rest  of  the  structure  in  the  form  of  stress  waves. 
Far  field  damage  away  from  the  Impact  area  could  result  from 
the  reflection  of  stress  waves.  It  is  generally  agreed  that 
the  cause  of  the  sudden  failure  must  be  examined  from  the 
point  of  transient  wave  propagation  phenomena. 

Flexural  waves  Induced  by  dynamic  loads  in  laminated 
composites  are  more  complicated  than  those  in  homogeneous 
and  isotropic  plates  due  to  the  anisotropic  and 
nonhomogeneous  properties  in  the  laminate.  Moreover, 
because  of  the  low  transverse  shear  modulus  In  fiber 
composites,  the  effect  of  transverse  shear  deformation 
becomes  significant  and  should  be  considered  in  the 
formulation.  In  Chapter  2,  the  laminate  theory  which 
Includes  the  transverse  shear  deformation  effect  is 
reviewed,  and  harmonic  waves  in  a  graphite/epoxy  laminated 
plate  are  studied.  The  propagation  of  wave  front  which,  for 
a  given  time  after  impact,  bound  the  stressed  region 
surrounding  the  impact  point,  is  also  Investigated. 

A  survey  of  wave  propagation  and  Impact  in  composite 
materials  has  been  given  by  Moon  [1].  Many  analytical  [2- 
5],  numerical  [6-7]  and  experimental  [8-10]  methods  have 
been  employed  to  study  the  transient  impact  problems.  The 
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respone  of  a  laminated  plate  can  be  analyzed  using  these 
methods  provided  the  applied  load  history  is  prescribed. 
Howfiver  if  the  dynamic  load  results  from  an  impact  of  an 
objcict  on  the  laminated  plate,  then  the  resulting  contact 
force  must  be  determined  first.  An  accurate  account  of  the 
concact  behavior  becomes  the  most  important  step  in 
analyzing  the  impact  response  problems. 

A  classical  contact  law  between  two  elastic  spheres  was 
derived  by  Hertz  [11].  When  letting  the  radius  of  one  of 
the  spheres  go  to  infinity,  one  obtains  the  contact  law 
between  an  elastic  sphere  and  an  elastic  half-space.  Many 
autliors  have  used  the  Hertzian  contact  law  for  the  study  of 
impact  on  metals  and  composites  [12-13].  Recently,  Yang  and 
Sun  [14]  performed  statical  indentation  tests  on  graphite/ 
epoxy  composite  laminates  using  spherical  steel  i ndenters  of 
dif-erent  sizes  and  found  that  the  Hertzian  law  of  contact 
was  not  adequate.  In  particular,  they  found  that 
sigtiif leant  permanent  indentations  existed  and  that  the 
unloading  paths  were  very  different  from  the  loading  path. 
Noting  that  energy  dissipation  takes  place  during  the 
process  of  impact,  Yang  and  Sun  [14]  suggested  that  this 
energy  is  responsible  for  the  local  damage  of  the  target 
materials.  The  unloading  curves  and  permanent  indentations 
obtained  from  the  statical  indentation  tests  may  provide  a 
use-'ul  information  in  estimating  the  amount  of  damage  due  to 
impact  since  this  energy  is  simply  the  area  enclosed  by  the 
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loading-unloading  curves.  In  this  study,  similar  statical 
indentation  tests  were  conducted  and  the  results  are 
presented  in  Chapter  3. 

Wang  [15]  has  performed  a  number  of  Impact  tests  on 
graphite/epoxy  laminated  beams  and  plates.  It  was  shown 
that  the  strain  responses  calculated  using  finite  element 
metliod  and  the  statically  determined  contact  laws  from  [14] 
agr<3ed  with  the  experimental  measurements  quite  well.  This 
indicates  that  the  statical  Indentation  law  Is  reasonably 
adequate  in  the  dynamical  Impact  analysis.  It  was  also 
suggested  that  the  contact  force  should  be  measured 
experimentally  to  provide  an  additional  basis  for  comparison 
with  the  finite  element  solution  which  could  allow  further 
evaluation  the  applicability  of  the  contact  laws  in  impact 
analysis.  Chapter  4  describes  an  impact  experiment  on 
graphite/epoxy  laminated  plate  using  an  impact-force 
transducer  with  a  spherical  steel  cap  as  the  Impactor.  The 
contact  force  history  and  strain  responses  at  various  points 
on  the  plate  were  measured  by  means  of  the  transducer  and 
surface  strain  gages,  respectively,  and  were  compared  with 
the  predictions  of  finite  element  analysis  using  the 
statically  determined  contact  law. 

Chapter  5  summarizes  the  results  obtained  in  Chapter  2,  3 


and  4. 
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CHAPTER  2 

STRESS  WAVE  IN  A  LAMINATED  PLATE 

A  laminated  plate  theory  which  includes  the  effects  of 
transverse  shear  deformation  and  rotatory  Inertia  was 
developed  by  Yang,  Norris  and  Stavsky  [16]  in  a  way 
suggested  by  Mindiin  [17]  for  homogeneous  isotropic  plates. 
It  was  shown  that  the  frequency  curves  for  the  propagation 
of  harmonic  waves  In  an  Infinite  two- layer  isotropic  plate 
In  plane  strain  agreed  with  the  predictions  of  the  exact 
solution  obtained  from  theory  of  elasticity  very  well.  A 
similar  laminated  plate  theory  was  developed  by  Whitney  and 
Pagano  [18]  and  was  employed  in  the  study  of  static  bending 
and  vibration  for  antisymmetric  angie-piy  composite  plates 
with  particular  layer  properties.  It  was  found  that  the 
effect  of  shear  deformation  can  be  quite  significant  for 
composite  plates  with  span-to-depth  ratio  as  high  as  20. 
Good  agreement  was  also  observed  In  numerical  results  for 
plate  bending  as  comparing  with  exact  solutions  of 
elasticity.  In  this  study,  the  laminate  theory  developed  by 
Whitney  and  Pagano  was  used  for  its  simplicity  yet  quite 
satisfactory  in  describing  the  harmonic  wave  propagation 
[19]. 


6 


2.1  Laminate  Theory  with  Transverse  Shear  Effects 

2.1.1  Lamina  Constitutive  Equations 

A  laminated  plate  of  constant  thickness  h  consists  of  a 
number  of  thin  lamlnas  of  unidirectional ly  f Iber-relnforced 
composite  perfectly  bonded  together.  Each  lamina,  whose 
fiber  may  orient  In  any  arbitrary  direction,  can  be  regarded 
as  a  homogeneous  orthotropic  solid.  Consider  a  typical  k— th 
lamina.  A  coordinate  system  (xi ,  Xa,  Xg)  is  chosen  In  such 
a  way  that  the  x^— Xa  plane  coincides  with  the  midplane  of 
lamina,  and  x,  and  Xa  axes  are  parallel  and  perpendicular  to 
the  fiber  direction,  respectively.  If  a  state  of  plane 
stress  parallel  to  the  x^—Xa  plane  is  assumed,  then  the  in— 
plane  stress-strain  relations  are  given  by 

k  k 

1 1  Qi  1  Qi  2  0  £  1 1 

^22  '  “  Ql 2  Q22  ®  '  ®22  ■  (2—1) 

“T 1  a-  ®  Qe  6.  .'^12. 

The  transverse  shear  stress— stra 1 n  relations  are  given  by 


(2-2) 


In  which 
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Ql1  — 

Q22  ~  E2  /  (  1  “V  1  2  2  1  ) 

Ql  2  ~  1  2^2  /  ^  ^  1  2 ^  2  1  ^  ~  ^  2  1  ^  1  /  ^  ^  1  2  ^  2  1  ) 

066=^12  (2-3) 

Q44  =  G23 
Qb  S  “  3 

are  the  so-called  reduced  stiffnesses,  where  E,  G  and  v  are 
Young's  modulus,  shear  modulus  and  Poisson's  ratio, 
respectively,  and  subscripts  1  and  2  denote  the  directions 
parallel  to  Xi  and  X2  axes,  respectively. 

The  coordinate  system  for  an  arbitrarily  oriented  lamina 
does  not,  in  general,  coincide  with  the  reference  axes 
(x,y,z)  of  laminated  plate  (see  Figure  2.1).  Using  the 
coordinate  transformation  laws  for  stress  and  strain,  we 
obtain  the  stress-strain  relations  in  laminate  reference 
system  as 

k  _  _  k 

C^xx  Q1I  Qi2^16®  ®  ^xx 

O’ y  y  Ql  2  ^22  Q26  Q  Q  ^yy 

■  txy  •  =  Ql  6  Q26  Qee  Q  Q  ■  '*^xy  ■  (2—4) 

Ty2  0  0  0  Q44  Q45  "^yz 

txz.  0  0  0  Q45  Qsb,  Yxz. 

in  which  ^ j  j  are  g i ven  by 

^11  —  Ql  1 01^+2 ( Ql  2'*'2Q6  6  )oi^n^+Q2  2^^^ 


9 


^22  ~  Qi  1  +2 ( Qi  2'*"2Q5 e ) n^+Q2 2*^^ 

^12  ~  ^Ql1  "^02  2"'4Q6  ^ 

^16  ~  (  Ql  1  ~Ql  2“2Q6  6  (  Ql  2"“Q2  2'^2Q6  6 

^26  ~  ^Ql1  ~Ql  2  ”2^6  6  ^  inn^+  (  Qi  2""Q2  6  )  fTl®  M  (  2“5  ) 

^66  ~  (Ql  1  ■I’Qa  2~2Qi  2~2Q0  6  6  ^ 

^44  ~  Q4  4*T>^+Q5  5 

^4  5  =  (Q4>4-Q5  5)mn 

^55  ~  Q4  4l^^'*'Q5 

where 

m  =  COS0  n  =  stn0 

and  0  !s  the  angle  between  x-axIs  and  x^-axls  measured  from 
X  to  Xi  counterclockwise  as  shown  In  Figure  2.1. 


2.1.2  Plate  Strain-Displacement  Relations 

The  displacement  components  of  the  laminated  plate  are 
assumed  to  be  of  the  form  [16] 

u(x,y,z)  =  u°(x,y)  +  z<P^ix,y) 

v(x,y,z)  =  v°(x,y)  +  z<^y(x,y)  (2-6) 

w(x,y,z)  =  w“(x,y)  =  w(x,y) 

where  u°,  v°  and  w°  are  the  midplane  displacement  components 
in  the  X-,  y-  and  z-d i rect i ons ,  respectively,  and  ‘^v 

are  rotations  of  cross-sections  perpendicular  to  x-  and  y- 
axis,  respectively  (see  Figure  2.2).  In  Equation  (2.6)  we 
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Figure  2.2  Laminate  displacement  components  for  a  cross- 
section  perpendicular  to  the  y-axis 


have  assumed  that  u  and  v  vary  linearly  In  the  thickness 
direction,  while  w  Is  constant  through  the  thickness. 

The  strain  components  for  a  point  In  k-th  lamina  of  the 
laminated  plate  with  a  distance  z  from  the  midplane  can  be 
computed  as 

^xx"'  =  ex°  + 

Cyv''  =  ey°  +  ZACy 

•V,,y'‘  =  V,y'’+  (2-7) 

Vyz''  =  3w/9y  +  9v/9z  =  9w/9y  +  <t>y  =  'Vy^® 

'Vxz'*  =  9w/9x  +  9u/9z  =  9w/9x  +  <j>y^  =  “Vxz” 

where 

•Vx®  =  9u°/9x 

■Vy®  =  9v°/9y  (2-8) 

ixy°=  9u°/9y  +  9v°/9x 

are  the  In-plane  strain  components  of  ml  dp  lane,  and 
ACx  =  d<P^/dx 

Ky  =  9<^y/9x  (2-9) 

ACxy^  9^x/9y  +  9^y/9x 

are  the  rotation  gradients. 

In  Equation  (2-7),  since  w,  </>x  and  <py  are  Independent  of 
z.  It  follows  that  the  transverse  shear  strains  are  constant 
through  the  thickness  of  the  plate. 
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Equation 

Thus,  the  strain  components  at  any  point  in  the  plate  can  be 
determined  from  the  extensional  strain  components  of  the 
midpiane,  the  rotation  gradients  of  the  plate  and  the 
distance  z  from  the  midpiane. 


(2-7)  can  be  written  in  concise  matrix  form  as 


k 


2.1.3  Stress-Resultants  and  Laminate  Constitutive  Equations 

Substitution  of  Equation  (2-10)  In  Equation  (2-4)  gives 
the  stress  components  for  a  point  in  the  k-th  lamina  as: 


The  stress-resultants  acting  on  a  laminate  can  be 
obtained  by  Integration  of  the  stresses  in  each  lamina 
through  the  laminate  thickness.  Specifically,  the  in-plane 


stress-resultants  are  given  by 


(2-12) 


the  stress  couples  are  given  by 


(2-13) 


and  the  transverse  shear  forces  are  given  by 


(2-14) 


The  sign  convention  for  these  stress-resultants  along 
with  the  geometry  of  a  typical  N-layer  laminated  plate  are 
shown  in  Figure  2.3. 

Substituting  Equation  (2-11)  into  the  right  hand  sides  of 
the  above  three  equations  and  performing  the  i ntegrat ions , 


we  obtain 
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where 

(Ai  j.Bij  ,Dij)  =  f j(1  ,z,z2)dz  i,J  =  1,2,6  (2-18) 

J-h/2 

and 

A*ij  =  f  ^ijdz  i,J  =  4,5  (2-19) 


Equations  (2-15)  through  (2-17)  are  usually  written 
symbol  leal ly  as 

A  B  0  1 [e° 

B  D  0  K  -  (2-20) 

0  0  A*J  lY 

which  Is  the  laminate  constitutive  equation  with  transverse 


shear  effect  Included. 
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2.1.4  Plate  Equations  of  Motion 

The  stress-equations  of  motion  for  the  k-th  lamina  are 
g 1 ven  by 


X  *  X 

4* 

X  V  J 

)  y 

+ 

^  X  z  >  z 

=  pu 

^  X  V  j  X 

+ 

Cf  y  Y  ! 

'  y 

+ 

^  y  z  >  z 

II 

<: 

X  Z  >  X 

+ 

^  y  z  J 

>  V 

+ 

^  z  z  >  z 

=  pw 

(2-21 ) 


where  p  Is  the  mass  density.  Integrating  Equation  (2  21) 
through  the  thickness  of  laminate  and  then  substituting 
Equation  (2-12),  (2-14)  and  (2-6)  in,  we  obtain 


Nx.x  +  =  Pu°  +  R^x 

Nxv»x  +  Ny,y  =  Pv°  + 

Qx  ,  X  +  Qy  .  y  +  q  =  Pw 


(2-22) 


where  q  is  the  normal  traction  on  the  plate.  Multiplying 
the  first  two  equations  of  Equation  (2-21),  integrating 
through  the  thickness  of  laminate  and  then  substituting 
Equations  (2-13),  (2-14)  and  (2-5)  in,  we  obtain 


Mx.x  +  M,y.y  -  Qx  =  Ru°  + 
Mxyix  «  V  Qy  “  RV°  + 

in  which  P,  R  and  I  are  defined  as 


(2-23) 


Jh  /  2 

p(1 ,z,z2)dz 

-h/  2 


(2-24) 


Equations 


(2-22)  and  (2-23)  are  the  plate  equations  of 
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motion.  Substitution  of  Equation  (2-20)  and  then  the 
strain-displacement  relations  In  these  two  equations  yield 
the  equations  of  motion  In  terms  of  midplane  displacements 
and  rotations  of  the  plate. 

A  graph  I te/epoxy  laminated  plate  provided  by  NASA  Lewis 
Research  Center  was  used  throughout  this  study.  This 
laminate  Is  a  [0°/45°/0°/-45®/0°] as  graph  I te/epoxy  composite 
with  0.0053  inch  ply  thickness  and  the  following  ply 
properties  [15]: 

El  =  17.5  X  10®  psi . 

Ea  =  1 . 15  X  10®  psi . 

Gia  =  Gi3  =  G23  =  0-8  X  10®  psi.  (2-25) 

Via  =  0.30 

p  =  0.000148  lb-sec^/ In'* 

For  symmetrically  laminated  composite  plate,  B|j  =0  and 
R  =0.  In  addition,  by  choosing  the  x-axis  of  the  laminate 
reference  system  to  coincide  with  the  0°  fiber  direction,  we 
obtain  Ai 6  =  Aae  =  0  and  Di g  =  Dae-  Further,  In  this  study, 
we  assume  0,3  =  Gas  =  G12.  and  consequently,  A* 45  =  0  and 
A* 4 4  =  A*55.  For  this  particular  laminate,  the 
displacement-equations  of  motion  are  given  by 

Ai,a2u®/3x=  +  AgeS^uVay^  +  (A,a  +  A66)9'^vV9x9y  =  Pu° 


(A, a  +  A66)9*u°/9xay  +  /dx^  +  Aa29^v°/ay^  =  Pv® 
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Dii92<#>x/8x2  +  2Di63^<^x/9x0y  +  D669^<^x/9y^ 

+  Di6^9^^y/9x^  +  9^^y/9y^)  +  (D^  +  De  s )  9^^y/9x9y 
-A*44(9w/9x  +  </>x)  =  (2-26) 

D,6(9=^x/9x2  +  9<#>x/9y^)  +  (Di2  +  D6  6)9^<#>x/9x9y 
+  De69  =  <#>v/9x2  +  2Di6  9^<#>v/9x9y  + 

“A*4  4(9w/9y  +  </>y)  =  I^y 

A*44(0=w/9x2  +  9=w/9y=  +  90x/9x  +  9<#>y/9y)  +  q  =  Pw 

In  Equation  (2-26),  the  first  two  equations  govern  the 
In-plane  motion  whi le  the  last  three  equations  govern  the 
flexural  motion. 

2.2  Propagation  of  Harmonic  Waves 

Consider  a  Infinitely  large  laminated  plate  governed  by 
the  equations  of  motion  (2-26).  We  assume  plane  harmonic 
waves  I n  the  form 

u°  =  U  exp[Ik(77  -  ct)] 
v°  =  V  exp[il<(r7  -  ct)] 

w  =  W  exp[ik(r7  -  ct)]  (2-27) 

=  ^x  exp[ik(77  -  ct)] 

^y  =  Sy  exp[lk(77  -  ct)] 

propagating  over  the  plate,  where  U,  V,  W,  ®x  and  5y  are 
constant  amplitudes,  k  is  the  wave  number,  c  i s  the  phase 
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velocity  and  77  Is  given  by 

77  =  X  cosa  +  y  sina  (2-28) 

in  which  a  is  the  angle  between  the  direction  of  wave 
propagation  and  x-axis. 

Substitution  of  Equation  (2-27)  into  Equation  (2-26)  with 
q  =  0  yields  a  system  of  five  homogeneous  equations  for  the 
five  constant  amplitudes.  In  order  to  have  a  nontrivial 
solution,  the  determinant  of  the  coefficient  matrix  is  set 
equal  to  zero.  Since  the  equations  are  uncoupled  into  two 
groups,  the  determi nantal  equation  can  be  seperated  into  two 
equations  as 

|a,jl=0  (2-29) 

for  the  in-plane  extenslonal  and  in-plane  shear  waves,  and 

|b,j|=0  (2-30) 

for  the  flexural  waves.  In  Equations  (2-29)  and  (2-30)  the 
coefficients  ajj  and  bjj  are  given  by 

ai 1  =  Aiicos^a  +  AesSln^a  -  Pc^ 

a-]2  —  ^21  —  (^12  Agg)si  nofcosoi  ( 2—3 1 ) 

£>22  =  Aggcos^a  +  A22Sln^cd  -  Pc^ 

and 

bi  1  =  Diik^cos^oc  +  2Di  gk^sinacosa  +  Dggk^sin^a 

+  A*44  -  Ik=c2 
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bia  “  bai  =  Diek^cos^a  +  (Dia  +  Dge )k^s I nacosa 
+  Di gk^sln^a 

bi3  =  bg,  =  iA*44kcosa  (2-32) 

baa  =  Deek^cos^a  +  20, ek^s i nacosa  +  Daak^sin^a 
+  A*44  -  Ik^c^ 

bas  =  bga  =  IA*44ksInof 

bga  =  -A*44k2  +  Pk^c^ 

Expanding  Equation  (2-29)  we  obtain  a  quadratic  equation 
in  as 

-  dic2  +  da  =  0  (2-33) 


where 


di 


da 


(AiiCos^a  +  Aaasin^a  + 


(2-34) 


Aiicos^a  +  Aeesin^of  (A,  a  +  Age )s i nofcosa 
(A,  a  +  A66)s'>^“<^Osa  AeeCOS^a  +  AaaSin^a 


It  is  noted  that  the  phase  velocity  c  does  not  depend  on 
the  wave  number  k,  thus  these  waves  are  nondi spers i ve.  In 
studying  of  transverse  impact  problem  where  in-plane 
deformation  is  negligible,  this  nondispersi ve  property  has 
no  significant  effect.  Should  In-plane  deformation  become 
important,  higher  order  approximation  of  displacement 
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components  u  and  v  must  be  assumed  and  the  dispersive 
property  of  these  in-plane  waves  could  be  included. 

From  Equation  (2-34)  It  is  evident  that  there  exist  two 
phase  velocities  corresponding  to  two  modes  of  wave. 
Although  these  two  waves  involve  both  in-plane  extensional 
deformation  as  well  as  in-plane  shear,  from  the  eigenvectors 
we  are  able  to  tell  which  one  is  dominant.  Thus  we  label 
the  two  waves  as  in-plane  extensional  wave  and  in-plane 
shear  wave  accordingly. 

The  determi nantal  equation  given  by  Equation  (2-30) 
yields  three  positive  roots  In  c^  Indicating  that  three 
flexural  waves  exist.  These  phase  velocities  are  functions 
of  the  wave  number  k,  thus  they  are  dispersive.  Among  these 
three  modes  of  wave,  only  the  lowest  one  corresponding  to 
the  transverse  shear  wave  has  a  finite  velocity  as  k—0  or  as 
the  wave  length  becomes  infinite.  The  dispersion  curves  for 
the  waves  of  the  lowest  mode  propagating  in  the  directions 
of  0°,  45°  and  90°  respectively  are  plotted  in  Figure  2.4 
with  the  non-dimensional  phase  velocity  vs.  the  non- 
dimensional  wavelength  X/h.  It  can  be  seen  that  they  all 
approach  the  value  of  yG^ 3/p  as  the  wavelength  becomes 
shorter.  The  phase  velocities  for  the  two  higher  modes, 
however,  approach  different  values  in  different  propagation 
directions  when  X-*0.  For  laminated  composite  which  are 
anisotropic  in  general,  the  phase  velocity  varies  from  one 
direction  to  another.  As  a  result  the  wave  surface  will 


C//  G12/  P 
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become  a  rather  complicated  shape  as  it  propagates.  This 
will  be  discussed  in  the  next  section. 

Substitution  of  o)  =  kc  In  Equation  (2—32)  yields  a  set  of 
frequency  equations  for  flexural  waves..  Figure  2.5  shows 
the  frequency  curves  of  these  waves  for  a  =  0° ,  45°  and  90°, 
respectively,  with  the  non-dimensional  frequency  vs.  the 
non-dimensional  wavelength.  The  cutoff  frequencies  for  the 
two  higher  modes  have  a  value  of  71 20, 3/p/h.  Comparing  with 
the  exact  cutoff  frequency  (7T/h)\/Gi  3/p ,  it  can  be  seen  that 
if  the  shear  correction  factor  7r^/12  is  introduced,  this 
theory  will  predict  the  correct  cutoff  frequency. 


2.3  Propagation  of  Wave  Front 

Impact  of  foreign  objects  on  a  laminated  plate  with  a 
very  short  duration  could  generate  weak  shock  waves  which 
will  propagate  Into  the  rest  of  the  structure  with  finite 
velocities,  and  the  positions  of  the  wave  fronts  define  the 
regions  being  disturbed  at  any  instant  after  impact. 
Damages  to  the  laminated  plate  may  possibly  occur  as  the 
first  wave  front  hits  the  weakest  part.  It  is  hence 
important  to  investigate  the  propagation  of  these  shocks  in 
the  plate.  There  have  been  works  dealing  with  the 
propagation  of  wave  front  in  anisotropic  elastic  media  [20— 
22].  Moon  [23]  presented  an  analysis  of  wave  surfaces  in  a 
laminate  by  treating  it  as  an  equivalent  homogeneous 
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orthotropic  plate.  The  acceleration  waves  and  their  wave 
fronts  were  Investigated.  The  propagation  of  shock  waves  In 
more  general  laminates  which  exhibit  the  bend Ing-ex tens lonal 
coupling  were  studied  by  Sun  [2].  The  ray  theory  was 
employed  to  construct  the  wave  front  surface.  The  growth 
and  decay  of  the  shock  strength  were  also  discussed.  In 
this  section,  the  analytical  procedures  developed  by  Sun  [2] 
were  applied  on  a  [0°/45VOV-45VOO]  as  graph  1  te/epoxy 
laminated  plate. 


2.3.1.  Kinematic  Conditions  of  Compatibility  on  the  Wave 
Front 

A  wave  front',  which  will  be  denoted  by  n,  Is  defined  as  a 
surface  travelling  over  the  plate  as  time  varies 
continuously,  and  across  which  there  may  exist  a 
discontinuity  In  the  stress,  particle  velocity  and  their 
der i vat  1 ves . 

Consider  a  discontinuous  surface  n  passing  some 
observation  point  In  a  medium  at  a  certain  time  t.  Let  f 
be  the  value  of  a  field  function  f(Xi ,t)  (e.g.  stress, 

particle  velocity,  etc.)  behind  the:  surface  fi,  and  f  be 
the  value  of  f  In  front  of  It,  then  the  discontinuity  of 
function  f  can  be  expressed  as 

[f]  = 


c  + 
\ 


f- 


(2-35) 
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In  which  the  right  hand  side  is  to  be  evaluated  at  the  time 
and  location  on  fl  passing  the  observation  point,  and  the 
Jump  across  the  wave  front  is  denoted  by  square  bracket. 

Surface  Q  may  be  expressed  mathematically  by  an  equation 


of  the  form 


5b(Xi  ,  t)  =  0 


(2-36) 


or,  by  making  t  explicit,  as 


5k(x,  ,t)  =  F(Xi)  -  t  =  0 


(2-37) 


which  represents  a  family  of  surfaces  in  Xj-space  with  t  as 
a  parameter.  By  evaluating  f^  and  f  at  t  =  F(x , ) ,  the  jump 
of  f  across  the  wave  front  becomes 


[f(x,)]  =  fMXi,F(Xi))  -  f-(Xi,F(Xi)) 


(2-38) 


The  rate  of  change  of  [f]  for  an  observer  moving  with  n 
Is  g I ven  by 

d[f]/dt  =  Of"/3x,  -  ef"/9Xi  )dXi /dt  +  {df*/dt  -  af/9t) 

=CiC9f/9Xi]  +  [9f/9t]  (2-39) 


v^here  t  =  F(X|)  is  substituted,  and  Cj  =  dxj/dt  are  velocity 
components  of  wave  front  relative  to  the  material. 

If  the  laminate  theory  introduced  in  previous  section  is 
used,  then  the  plate  displacement  components  are  u° ,  v°,  w, 
and  4>y,  while  the  spatial  variables  are  Xi  =  x  and  Xa 
y.  It  is  assumed  that  the  Integrity  of  the  material  is  not 
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affected  by  the  propagation  of  the  stress  wave  front,  then 
these  displacement  components  will  remain  continuous. 
Consequently,  we  have 

[u°]  =  [v°]  =  [w]  =  [<^J  =  i<f>yl  =  0  (2-40) 


across  the  wave  front.  Applying  the  general  condition  of 
Equation  (2-39)  on  u°,  together  with  Equation  (2-40),  we 
obtain 


[auV9Xj]cj  +  [u°]  =0  j  =  1  ,2 


(2-41) 


Let  Cn  and  nj  be  the  normal  velocity  and 
on  the  wave  front,  respectively,  it  follows 


the  un i t 
that 


norma  1 


HjCj  =  Cn 

and  Equation  (2-41)  becomes 

[au°/axj]  =  -[u“]nj/Cn  J  =  1,2 


(2-42) 


(2-43) 


Similar  relations  can  be  derived  for  the  other 
displacement  components  v° ,  w,  and  .  Together  they 
specify  the  kinematic  conditions  of  compatibility  on  the 
wave  front. 


2.3.2  Dynamical  Conditions  on  the  Wave  Front 

Consider  a  finite  volume  V  of  a  material  medium  and 
denoted  by  S  the  boundary  or  surface  of  V.  For  a  continuous 
and  differentiable  function  f(Xi,t)  inV,  it  can  be  shown 
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[23]  that 


,J  f(Xi,t)dV  =  +  J^GfdS 


dt J  V 


(2-44) 


under  deformation  of  the  medium,  where  G  is  the  normal 
veiocity  of  the  surface  S.  In  the  case  where  the 
deformation  of  the  volume  V  is  determined  solely  by  the 
motion  of  material  particles,  we  have 


G  =  Uini  =  Un 


(2-45) 


where  u,  is  the  displacement  components;,  n,  is  the  outward 
normal  on  S,  and  Un  is  the  normal  velocity  of  material 
particle  on  S.  If  there  exists  a  discontinuity  surface  (or 
wave  front)  travelling  with  velocity  Cj  in  the  medium,  by 
choosing  this  surface  as  the  boundary  of  V,  we  have 


G  —  c  j  n  j  Cn 


where  Cn  is  the  normal  velocity  of  wave  front. 


(2-46) 


Suppose  that  a  volume  V  whose  motion  is  determined  by  the 
deformation  of  the  material  medium,  is  divided  by  a 
travelling  surface  Q  into  two  volumes  V  and  y*  as  shown  in 
Figure  2,6,  The  surface  S  is  also  divided  into  two  portions 
S“  and  S'*"  which  form  parts  of  the  boundaries  of  V  and  V*, 
respectively.  The  remaining  part  of  the  boundary  is  formed 
by  Qo  which  Is  a  segment  of  Q.  In  Figure  2.6,  nj  denotes 
the  unit  normal  of  Q  in  the  direction  of  travelling,  and  n,* 
denotes  the  unit  outward  normal  of  S. 


30 


Taking  f  =  pu,  in  Equation  (2-44)  and  using  equation  (2- 
45)  and  (2-46),  we  obtain 


^J^puldV  =  +  J^^u^putdS  -  J^Cnputdn  (2 


CnpUidfi 

(2-47) 

Cnputdn 

(2-48) 

where  u7  and  ut  are  the  velocity  components  of  material 
particle  in  V  and  M* ,  respectively.  Combining  the  above 
two  equations  gives 

4|^pu,dV  .  J'^(pu,).tdV  +  J^0;puTdS  +  X^«pdtdS 


+  (*  Cnp(uT  -  u|)' 

V  fto 


(2-49) 


From  theory  of  elasticity  we  have 

-  ;^PunjdS  (2-50) 

If  we  let  the  volume  V  approach  zero  at  a  fixed  time  in  such 
a  way  that  it  will  pass  into  fio .  then  the  volume  integral  in 
Equation  (2-49)  will  evidently  approach  zero;  however 


J  «pC,Tds  -  -  J^^puT. 


(2-51  ) 
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f  u;pu7dS  -  r  UnpuTdn  (2-52) 

t/  s-  v  fto 

J^ffjjnjdS  -  “  <^Tj>njdn  (2-53) 

where  cr7j  and  atj  are  the  stress  components  on  the  sides  of 
Qo ,  respectively. 

Substituting  Equations  (2-50)  through  (2-53)  into 
Equation  (2-49)  gives 

J  (<^t  j“0'7  j  )njdn  =  J*  pu7(Cn-u")dn  -  J*^pu7  (Cn-Un ) dn  (2-54) 

Using  [ffij]  and  [Uj]  to  represent  the  Jumps  of  stress  and 
particle  velocity  across  the  wave  front,  and  utilizing  the 
fact  that  Cn  »  Un ,  we  obtain 


f  [ffijlnjdn  =  -  f  pCn[Ui]dn 

V  fto  ^ 


(2-55) 


Since  this  condition  is  independent  of  the  extent  of  the 
surface  integration  Qq ,  it  follows  that 

[(Tijlnj  =  -  pCnCu,]  (2-56) 

In  the  case  of  laminated  plate,  i  =  x,y,z  and  J  =  x,y. 

Substitution  of  Equation  (2-6)  into  Equation  (2-56) 


yields 
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[(Tijlrij  =  -  pCniCu®]  +  z[^>x3} 

[(Tajlnj  =  -  pc„{Cv°]  +  z[^v]} 

C<^3j]rij  “  “  pCnCw] 


(2-57) 


Integrating  Equation  (2-S7)  over  the  thickness  of  plate 


g  i  ves 


CNx]nx  +  ENxylny  =  -  Pcn[u°]  -  RcnC^x] 
[Nxy]nx  +  [NylOy  =  -  PCnCv°]  -  RCnC^y] 
CQxlnx  +  EQylny  =  -  PcnEw] 


(2-58) 


Multiplying  the  first  two  equations  of  Equation  (2  57)  by  z 
and  then  integrating  over  the  thickness,  we  obtain  two  more 


equat i ons 


[Mxlnx  +  [Mxy]ny  =  -  RCnCu°]  -  ICnC^xl 

CMxy]nx  +  CMy]ny  =  -  RcnCv°]  -  ICnE^yl 


(2-59) 


where  P,  R  and  I  have  been  defined  in  Equation  (2  24) 

The  five  equations  given  by  Equations  (2-58)  and  (2-59) 
are  the  dynamical  conditions  on  the  wave  front  for  the 

laminated  plate. 


2.3.3  Propagation  Velocity  of  the  Wave  Front 

Across  the  wave  front,  the  laminate  constitutive 
relations  given  by  Equation  (2-20)  can  be  written  as 
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CN] 
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where 


{[N]}^  =  {[Nj,[Ny].[N,y]} 
{[Q]}"^  =  {[Qxl.CQy]} 


(2-60) 


(2-61) 


are  the  Jumps  of  the  stress  resultants,  and 

{[e]}T  =  {  C9uVax]  ,  [3vV3y] ,  Cau°/ay]  +  [3v®/3x]  } 

{[k]}"^  =  {[a<#>x/3x],[30y/3y],[a<^x/9y]  +  C9<^x/9x]}  (2-62) 

=  { [3w/8y] ,  [aw/ax] } 

are  the  Jumps  of  the  strain  components.  In  Equation  (2—62), 
the  conditions  [<#>x]  =  =  0  substituted. 

Substituting  of  Equation  (2-43)  and  the  similar  relations 
for  other  kinematic  variables  in  Equation  (2-60),  we  can 
express  the  Jumps  of  the  stress  resultants  In  terms  of  the 
Jumps  of  the  time  derivatives  of  the  displacement  components 
u°,  v°,  w,  <3!>x  and  4>y .  These  relations  are  then  substituted 
in  Equations  (2-58)  and  (2-59),  which  results  in  five 
homogeneous  equations.  For  [0°/45°/0°/-45°/0°] 2s  graphite/ 
epoxy  laminated  plate  which  is  symmetrical  and  balanced 
(l.e.  Bij  =  0,  Aie  =  Aac  =0,  R  =  0  and  Die  =  Dae),  these 
five  equat I ons  are  uncoupled  Into  three  groups  as 
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(2-63) 

(2-64) 

(2-65) 


In  which  [a,jl  and  [b,j]  are  both  two  by  two  syrmetric 
matrices,  and  their  entries  are  given  by 


ai 1  =  nx^Aii  +  ny^Ass  " 

ai2  =  ^21  =  nxny(Ai2  +  Age) 

62  2  =  l^x^Age  +  ny^A2  2  “ 

bii  =  n,,2Dii  +  2n,,nyDi6  +  Hy^Dgs  “ 
bi2  =  ^21  =  Di6  +  nxny(Di2  +  Dge) 
b22  =  n^^^Dgg  +  2n,,nyD,6  +  ny=D22  “ 


(2-66) 


(2-67) 


It  can  be  seen  that  Equation  (2-63)  describes  the  in 
plane  extensional  and  the  in-plane  shear  wave  fronts, 
Equation  (2-64)  describes  the  bending  moment  and  the 
twisting  moment  wave  fronts  and  Equation  (2-65)  describes 
the  transverse  shear  wave  front. 

From  Equation  (2-65),  we  obtain  the  normal  velocity  with 
which  the  transverse  shear  wave  front  propagates  as 


2  _  A* 


A*44/P 


(2-68) 


It  is  noted  that  this  velocity  Is  Independent  of  the 

and  is  called  directionally 


direction  of  propagation, 
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isotropic  wave  front. 

Equations  (2-63)  and  (2-64)  yield  non-trivial  solutions 
only  if  the  determinant  of  the  coefficients  matrices  vanish, 

i  .e. 


la.jl  =  0 

b,j|  «  0 


(2-69) 

(2-70) 


Each  of  the  above  equations  can  be  expanded  into  a 
quadratic  equation  of  Cn^.  For  [0°/45°/0°/-45®/0°] as 
graphite/epoxy  laminated  plate,  the  normal  velocities  of 
wave  fronts  corresponding  to  the  in-plane  modes  and  flexural 
modes  are  plotted  in  Figure  2.7  and  2.8,  respectively.  It 
is  noted  that  the  normal  velocities  of  the  in-plane 
extenslonal  and  in-plane  shear  modes  are  symmetrical  about 
x-axis  and  y-axis,  while  there  is  no  such  symmetry  for  the 
bending  moment  and  twisting  moment  modes. 


2.3.4  Wave  Surface  and  Ray 

From  Figure  2.7  and  2.8,  it  can  be  seen  that  for 
laminated  composites  which  are  anisotropic  in  general,  the 
in-plane  and  flexural  wave  fronts  travel  with  different 
normal  velocities  in  different  directions.  In  other  words, 
the  initial  shape  of  a  wave  surface  will  be  distorted  after 
it  propagates.  However,  Equations  (2-66)  and  (2-67)  show 
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Figure  2.7  Normal  velocities  of  in-plane  wave  fronts 
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that  for  any  fixed  normal  direction  n, ,  Cn  Is  a  constant. 
Connecting  the  points  having  the  same  unit  normals  to  the 
travelling  wave  front  surface,  we  obtain  a  family  of  lines 
which  are  cal  led  rays.  Thus,  along  a  ray,  the  normal 
velocity  of  wave  front  remains  unchanged.  By  using  the  ray 
theory  which  has  been  well  established  In  the  field  of 
geometrical  optics,  we  are  able  to  construct  the  wave  front 

surface. 


Recall  Equation  (2.37) 


F(x, )  -  t  =  0  1=1,2 


(2-37) 


which  represents  a  family  of  wave  fronts  propagating  over 
the  plate  with  t  as  a  parameter.  It  follows  that 


dF/dt  =  (9F/aXi ) (dx,/dt)  =  (aF/aXi)ci  -  1 


(2-71 ) 


By  putting 


=  aF/ax,  =  VF 


(2-72) 


Equation  (2—71)  becomes 


PiCi  =  1 


(2-73) 


Since  Pi  Is  normal  to  the  surface  F,  It  can  be  written  as 


Pi  =  |Pi I 


(2-74) 


y/here  [Pi  |  denotes  the  length  of  Pi  . 
(2-74),  we  obtain 


Combining  (2-73)  and 
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|Pi |niC,  =  Ip, |Cn  =  1  (2  75) 

from  which  we  obtain 

Pi  =  Hi/Cn  (2-76) 

In  Equation  (2-76),  p,  is  called  the  slowness  vector 
which  has  the  direction  normal  to  the  wave  front  with  the 
magnitude  being  equal  to  the  inverse  of'  normal  velocity  Cn • 

Subsitituting  Equation  (2-76)  in  Equation  (2-69)  and  (2- 
70),  we  obtain  two  equations  in  terms  of  p, 

Px^Aii  +  Py^Ass  ”  P  PxPy(Ai2  +  Aec) 

PxPv^Aia  +  Ass)  Px^Aes  +  Py^Aaa  “  P 

p,,2Di  i+2p,,PyD,  e+Py^Dee-I  Di  6+PxPy(Di  a+Dse) 

Di 6+PxPy(Di s+Dee)  Px 6+2PxPyDi e+Py^Daa-I 

which  can  be  written  in  a  general  form  as 

g(pi)=0  i=1,2  (2-77) 

In  view  of  Equation  (2-72),  we  recognize  that  Equation 
(2-77)  may  be  regarded  as  a  set  of  first-order  partial 
differential  equation  for  F.  A  standard  method  of  solving 
first-order  partial  differential  equation  is  by  means  of 
characteristics  [24],  which  reduces  the  equation  to  a  system 
of  first-order  ordinary  differential  equations.  In  our 
case.  Equation  (2-77)  then  is  equivalent  to  the  following 


dx/ds  =  8g/apx  dy/ds  =  ag/8py 
dpx/ds  =  -9g/9x  dpy/ds  =  -9g/9y 
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(2-78) 
(2-79) 

vj,here  s  Is  a  parameter.  These  equat Ions , together  with 
Equation  (2-77)  describe  the  ray  geometry  and  the  normal 
direction  of  the  wave  front  propagating  along  the  ray. 


From  Equation  (2-78),  we  have 

dy/dx  =  (9g/9py)/(9g/9px)  ^2-80) 

Since  the  normal  direction  of  wave  front  along  a  ray  is 
constant,  It  can  be  seen  from  Equation  (2-76)  that  Pi  is 
also  constant  along  a  ray.  For  laminated  composite  which  is 
assumed  to  have  homogeneous  material  properties.  Equation 
(2-77)  shows  that  g(pi)  does  not  depend  on  Xj ,  consequently, 
ag/3px  and  9g/9py  are  all  constants  along  a  ray.  Thus,  the 
solution  of  Equation  (2—80)  is  then  given  by 

y  =  f (x  -  Xo)  +  yo  (2-81 ) 

where  Xo  and  yo  are  the  initial  values  of  x  and  y  at  t  =  0, 
and  f  =  (9g/9py)/(9g/9px).  This  equation  shows  that  the 
rays  In  a  homogeneous  solid  are  straight  lines. 


From  Equations  (2-73)  and  (2  77),  we  have 

Cjdpi  =0 


(2-82) 


dg  =  (9g/9pi)  dpi  =  0 


(2-83) 
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Eliminating  dpj  from  these  equations  yields 

dXi/dt  =  C|  =  (0g/dpi )/(Pj9g/9pj )  (2-84) 

where  summation  over  J  Is  understood. 

Equation  (2-84)  can  be  solved  to  obtain  the  position  of 
wave  front  at  time  t.  Again,  since  9g/9pj  and  Pj  are  all 
constant  along  a  ray,  we  obtain  the  solution  of  Equation  (2- 
84)  as 

X  =  (ag/9px)t/(pj9g/9pj )  +  Xq  (2-85) 

y  =  (ag/9py )t/(pj9g/9pj )  +  yo  (2-86) 

where  Xq  and  yo  denote  the  initial  wave  position  at  t  =  0. 

Consider  at  t  =  0,  a  wave  front  forms  a  circle  given  by 


Xq  =  h  COSOf 

yo  =  h  sina 


(2-87) 


At  this  instant,  the  normal  directions  to  the  wave  front 
coincide  with  the  radial  directions.  Due  to  the  different 
velocities  of  propagation  in  directions,  this  initial  shape 
would  be  distorted.  By  using  Equations  (2—85)  and  (2—86), 
the  subsequent  positions  of  the  wave  front  can  be 
determined.  Figures  2.9-2.12  show  the  wave  front  positions 
at  two  consecutive  instants  after  t  =  0  for  the  in-plane 
extensional,  In-plane  shear,  bending  moment  and  twisting 
moment  modes,  respectively,  for  the  [0°/45°/0°/-45°/0°] as 
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graph  I te/epoxy  laminated  plate.  It  is  noted  that  for 
symmetrical  laminates,  the  in-plane  modes  are  uncoupled  from 
the  bending  modes.  The  rays  along  which  the  normal 
directions  to  the  wave  front  are  0°.  45°  and  90°, 

respectively,  are  also  shown  in  the  figures.  It  is  seen 
that  the  wave  fronts  of  the  in-plane  extensionai  and  the  in¬ 
plane  shear  modes  possess  symmetry  with  respect  to  x~axis 
and  y-axis.  The  wave  fronts  of  the  bending  and  twisting 
moments,  however,  lose  their  original  symmetry  with  respect 
to  x-axis  and  y-axis.  This  is  an  indication  that  in 
performing  analysis  of  flexural  deformation  of  this 

laminate,  one  can  not  take  a  quadrant  for  analysis,  a 

practice  followed  by  many  authors  dealing  with  homogeneous 

and  isotropic  plates. 

From  Figures  2.9-2.12,  it  is  also  interesting  to  note 
that  ray  geometries  for  these  two  groups  of  wave  fronts  are 
quite  different.  For  the  In-plane  extensionai  and  in-plane 
shear  wave  fronts,  the  rays  coincide  with  the  normal 
directions  when  «  =  0°  and  90°.  Along  other  directions,  the 
direction  of  the  ray  deviates  from  the  normal  direction  of 
the  wave  front.  It  was  discussed  in  [2]  that  the  degree  of 
spreading  of  rays  is  proportional  to  the  decay  of  the  stress 
amplitude  at  the  wave  front.  Thus,  from  Figures  2.9  and 
2.11,  one  can  conclude  that  the  strength  of  the  in-plane 
extensionai  and  bending  moment  wave  fronts  decay  more 
rapidly  in  the  y-di recti  on  than  in  the  x-di recti  on. 
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Figure  2.9  Wave  front  positions  at  different  times  and 
rays  for  in-plane  extensional  mode 


Figure  2.10  Wave  front  positions  at  different  times  a 
rays  for  in-pianr  shear  mode 
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Figure  2.11  Wave  front  positions  at  different  times  and; 
rays  for  bending  mode 
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A  photoelastic  study  of  anisotropic  waves  In  a  fiber 
reinforced  composite  has  been  done  by  Dally  ^  [9]. 
The  waves  was  produced  by  a  explosive  charge  In  a  small  hole 
on  the  plate.  The  result  showed  clearly  an  el  1  I pt Ic-1 I ke 
stress  wave  front  pattern.  This  indicates  that  stress  waves 
in  anisotropic  materials  propagate  with  different  velocities 
in  different  directions. 
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CHAPTER  3 

STATICAL  INDENTATION  LAWS 

A  brief  introduction  of  the  historical  development  on 
impact  problem  involving  homogeneous  isotropic  materials  was 
given  by  Goldsmith  [12].  Hertz  [11]  was  the  first  to  obtain 
a  satisfactory  solution  on  contact  law  for  two  isotropic 
elastic  spherical  bodies.  When  letting  the  radius  of  one  of 
the  spheres  go  to  infinity,  this  law  then  describes  the 
contact  behavior  between  a  sphere  and  an  elastic  half-space. 
The  Hertzian  law,  In  spite  of  being  static  and  elastic  in 
nature,  has  been  widely  applied  to  impact  analyses  where 
permanent  deformations  were  produced.  The  use  of  this  law 
beyond  the  elastic  limit  has  been  Justified  on  the  basis 
that  it  appears  to  predict  accurately  most  of  the  impact 
parameters  that  can  be  experimentally  verified. 

In  studying  impact  responses  of  laminated  composites,  the 
problem  becomes  extremely  complicated.  One  may  easily 
realize  that  the  Hertzian  contact  law  which  was  derived 
based  on  homogeneous  isotropic  materials  may  not  be  adequate 
In  describing  the  contact  behavior  of  laminated  composites 
due  to  their  anisotropic  and  nonhomogeneous  properties. 
Moreover,  most  of  the  laminated  composites  have  finite 
thickness  which  can  not  be  represented  by  a  half-space.  In 
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many  existing  analytical  works  [25],  loadings  to  the 
laminates  were  assumed  known,  and  the  responses  of  the 
laminates  were  assumed  elastic. 

Willis  [26]  obtained  explicit  formulas  for  Hertzian 
contact  law  for  transversely  isotropic  half-space  pressed  by 
a  rigid  sphere,  and  extended  it  to  the  application  of  impact 
problems.  It  was  shown  that 

F  =  ka"  (3-1 ) 

with  n  =  3/2  is  valid  for  the  contact  force  F  and  the 
indentation  a,  where  k  is  a  contact  coefficient  whose  value 
depends  on  the  material  properties  of  the  target  and  the 
sphere,  and  the  radius  of  sphere. 

A  modified  contact  law  with 

k  =  (4/3)  - ^ -  (3-2) 

1  -  i/s""  1 

-  +  — 

Es  Et 

was  used  [13]  in  an  analytical  study  on  impact  of  laminated 
composites.  In  Equation  (3-2),  Rs.^s  Es  are  the  radius, 

Poisson's  ratio  and  Young's  modulus  of  the  sphere, 
respectively,  and  Et  is  the  Young's  modulus  of  the  laminates 
in  thickness  direction.  It  was  also  suggested  by  Sun  ^  al . 
[27]  that  the  value  of  k  can  be  experimentally  determined. 
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Recently  Yang  and  Sun  [14]  have  conducted  static 
indentation  tests  on  the  [00/45V0V-45V0°]  as  graphite/ 
epoxy  laminates  using  spherical  steel  indenters  of  0.25  in. 
and  0.5  In.  diameters.  The  results  were  fitted  into 
Equation  (3-1)  and  were  found  that  the  3/2  power  is  valid. 
In  addition,  it  was  also  observed  that  even  for  small 
amounts  of  load  there  were  significant  permanent 
Indentations.  This  Implies  that  the  unloading  curve  has  to 
be  different  from  the  loading  curves.  In  order  to  account 
for  the  permanent  deformation,  the  equation 

P  .  (3-3) 

\0!m  -  Wo  / 

proposed  by  Crook  [28]  was  used  to  model  the  unloading  path 
where  F„,  is  the  contact  force  at  which  unloading  begins,  a,, 
is  the  indentation  corresponding  to  F„,,  and  oio  denotes  the 
permanent  indentation  In  an  unloading  cycle.  Equation  (3-3) 
can  be  rewritten  as 


in  which 

s  =  F,„/(am  -  ao)*’ 

is  called  unloading  rigidity.  In  order  to  simpl ify  the 
modeling  of  the  unloading  law,  it  was  assumed  [14]  that  the 
for  all  the  unloading  curves  remains  the  same. 


value  of  s 
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Consequently,  a  constant  Wcr  given  by 
ttcr  =  k/s 


(3-6) 


was  introduced.  It  was  also  shown  that  q=5/2  fitted  the 
unloading  path  very  well,  and  the  permanent  Indentation  Wq 
was  then  related  to  Wm  by 


ao/“m  =  1  “  ttm  >  r 

=:  0  as  Dim  —  r 


(3-7) 


The  value  of  otcr  was  found  to  be  independent  of  the  size 
of  the  indenter  and  hence  can  be  regarded  as  a  material 

constant . 

It  was  also  mentioned  in  [14]  and  [29]  that  there  were 
some  practical  difficulties  in  performing  the  tests.  Since 
the  indentation  was  measured  step  by  step  using  a  dial  gage 
and  readings  on  the  gage  were  taken  about  10  to  20  seconds 
after  the  load  was  increased  by  one  step,  the  creep  effect 
may  cause  an  appreciable  error  to  the  results.  Another 
important  problem  was  that  it  was  almost  impossible  to 
measure  the  permanent  indentation  accurately  using  the  dial 
gage.  In  order  to  overcome  these  problems,  a  Linear 
Variable  Differential  Transformer  (LVDT)  was  used  in  this 
study  to  measure  the  indentation. 

The  LVDT  is  an  electromechanical  transducer  that  produces 
an  electrical  output  proportional  to  the  displacement. 
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Connecting  this  output  and  the  one  from  the  strain  indicator 
which  -is  used  to  measure  the  applied  ioading  to  a  X-Y 
:pl otter,  one  can  obtain  a  continuous  loading-unloading 
-curve.  By  changing  the  loading  rate  which  can  be  applied  as 
fast  as  50  ib./sec.,  it  is  possible  to  examine  the 
Significance  of  creep  effect  on  the  contact  law.  The 
starting  point  and  final  point  of  a  ioading-unloading  cycle, 
which  represent  respectively  the  instants  of  contact  and 
separation  of  the  indenter  and  the  specimen,  can  be  easily 
determined  from  the  curve.  Thus,  the  measurements  of 
permanent  indentations  are  much  more  accurate  than  those 
using  the  dial  gage. 


3.1  Specimens  and  Experimental  Procedure 

Two  groups  of  test  specimens  were  prepared  from  a  [0°/ 
450/00/-45VO°]2s  graphite/epoxy  laminate.  They  were  cut  In 
the  way  such  that  the  longitudinal  axis  of  the  beam  specimen 
of  the  first  group  was  para ii el  to  the  0°  fiber  direction 
while  the  second  one  was  perpendicular  to  it.  The  latter 
then  becomes  [90V45°/90V-45V90°]  laminated  beams.  The 
thickness  of  the  beam  was  0.106  In.  and  the  width  was 
approximately  1 .25  in. .  In  ail  tests,  the  specimens  were 
clamped  at  both  ends.  It  was  shown  in  [14]  that  the  span  of 
the  specimen  In  the  range  of  2  in.  to  6  In.  has  little 
effect  on  the  contact  law.  Hence,  only  one  span,  i.e.  2 
in.,  was  used  in  the  test. 


53 


The  experimental  set-up  Is  shown  schematically  In  Figure 
3.1.  LVDT  was  mounted  on  a  'C  bracket  fixed  to  the  loading 
piston  so  that  only  the  relative  movement  between  the 
indenter  and  the  specimen  was  recorded.  The  load  was 
applied  pneumaticallt  by  a  plunger  and  it  was  measured  using 
a  load  cell  and  a  strain  indicator.  Outputs  from  LVDT  and 
strain  indicator  were  fed  into  an  X-Y  plotter  so  that  a 
continuous  force- Indentation  curve  can  be  obtained.  Two 
spherical  steel  indenter s  of  diameters  0.5  in.  and  0.75  in. 
were  used . 


3.2  Experimental  Results 
3.2.1  Loading  Curves 

The  experimental  curves  were  first  digitized  into  some 
discrete  data  points  and  then  fitted  into  Equation  (3-1) 
using  least-squares  method.  Figures  3.2  and  3.3  show  the 
test  data  and  the  fitted  curves  for  0.5  in.  diameter 
indenter.  It  can  be  seen  from  these  figures  that  the  3/2 
power  index  gives  very  good  results.  However,  the  contact 
coefficient  k  of  [0V45V0°/-45°/0°]  as  specimen  is  less  than 
the  one  of  [90 V45V90°/-45 V90° ]  a s  specimen  by  about  7  %. 
During  the  test,  larger  deflections  were  observed  for  the 
second  group  of  specimen  due  to  their  lower  flexural 
rigidity.  This  means  that  the  contact  area  is  also  larger 
and  the  indentation  under  same  amount  of  loading  should  be 
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Figure  3.1  Schematical  diagram  for  the  indentation  test 
set-up 


K=1 .37S6E+06 
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smaller  comparing  with  the  first  group  of  specimens. 
Consequently,  the  higher  value  of  k  for  the  [90°/45°/90°/— 
45°/90°]2s  specimens  Is  reasonable. 

The  results  for  0.75  In.  diameter  Indenter  are  presented 
in  Figures  3.4  and  3.5.  Again,  good  agreement  between  the 
experimental  data  and  fitted  curves  indicates  that  the  3/2 
power  index  for  loading  law  is  val id.  The  values  of  k  for 
both  indenters  are  summarized  in  Table  3.1.  It  should  be 
noted  that  the  average  value  of  k  obtained  from  the  two 
groups  of  specimens  was  used  later  in  a  finite  element 
analysis  of  Impact  responses. 

3.2.2  Unloading  Curves 

By  choosing  a  suitable  value  for  q,  it  can  be  seen  from 
Equation  (3-5)  that  once  the  relation  between  oLq  and  is 
established,  the  unloading  rigidity  s  Is  then  determined. 
Test  results  show  that  the  permanent  indentations  ocq  and  the 
corresponding  maximum  indentations  (Xm  exhibit  a  rather 
linear  relationship.  The  equation  given  by 

=  Sp  (ttm  -  ttp)  (3-8) 

is  obtained  from  the  test  data  for  both  0.5  in.  and  0.75 
In.  indenters  using  least-squares  fitting  method,  and  are 
plotted  in  Figure  3.6.  In  Equation  (3—8),  (Xp  can  be 
considered  as  a  critical  value  of  indentation.  Once  the 
amount  of  Indentation  exceeds  (Xp,  permanent  deformation  will 
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These  two  equations  along  with  Equation  (3-4)  are  then  used 
to  fit  the  experimental  unloading  curves  in  finding  the 

value  of  q. 

Yang  [14]  has  shown  that  q  =2.5  fits  the  test  results 
for  both  0.25  In.  and  0.5  In.  i ndenters  qui te  wel 1 .  In 
this  study,  however,  the  values  of  2.2  and  1.8  were  found  to 
give  the  best  fitting  for  0.5  in.  and  0.75  in.  indenters. 
respectively  using  the  aforementioned  method  (Figures  3.7- 
3.10).  For  convenience,  q  =  2.5  was  used  for  0.5  in. 
indenter  whl 1e  q  =  2.0  was  chosen  for  3/4  in.  indenter. 
The  results  of  the  curve-fitting  are  presented  in  Figures 
3.11-3.14.  Further  discussions  on  the  unloading  law  will  be 
in  Section  3.3. 
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3.2.3  Reloading  Curves 
The  equation 

F  =  k,  (a  -  ao)P  (3-11) 

suggested  by  Yang  [14]  was  used  to  model  the  reloading 
curve,  where  k,  is  called  reloading  rigidity  and  p  =  3/2  was 
found  to  fit  the  experimental  data  quite  well.  It  was  also 
observed  that  the  reloading  curve  always  returns  to  where 
the  unloading  began,  and  hence  the  reloading  rigidity  can  be 
determined  by 

k,  =  F„,/(a,n  -  (3-12) 

In  other  words,  the  reloading  test  is  not  necessary  provided 
the  unioading  condition  is  specified.  Some  reloading  curves 
obtained  following  Equations  (3-11)  and  (3-12),  and  the 
experimental  data  are  presented  in  Figures  3.15-3.18. 


3.3  Discussion 

As  mentioned  before,  due  to  creep  the  loading  rate  may 
affect  the  contact  law  (l.e.  the  value  of  k) .  A  series  of 
tests  with  different  loading  rates  was  performed  to  examine 
this  point.  The  maximum  loading  rate  the  test  equipment  can 
apply  without  exceeding  it's  capacity  is  about  50  Ib/sec.. 
It  was  found  that  in  the  range  of  5  Ib/sec.  to  50  Ib/sec., 
the  values  of  k  showed  very  little  scatter,  and  the  effect 
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due  to  local  material  nonhomogeneity  in  the  composite  may  be 
even  greater  than  the  one  due  to  the  loading  rate.  However, 
an  appreciable  decrease  of  the  value  k  was  observed  when  the 
loading  rate  was  lower  than  1  Ib/sec..  In  some  extreme 
cases  where  loadings  were  applied  as  slow  as  10  Ib/min.,  the 
average  value  of  k  for  0.5  in.  i ndenter  was  very  close  to 
the  one  obtained  previously  by  Yang  [14]  using  dial  gage  to 
measure  the  ’indentation.  In  this  study,  the  loading  rates 
for  all  tests  were  approximately  equal  to  10  Ib/sec.. 

Unlike  the  exponent  n  of  the  loading  law  for  which  the 
value  of  3/2  seems  to  yield  good  agreement  with  all 
experimental  data,  the  exponent  q  of  the  unloading  law 
(Equation  3-3  or  3-4)  reveals  much  wider  deviation  for 
different  sizes  of  i ndenter.  Value  of  q  =  3/2  corresponding 
to  an  elastic  recovery  according  to  the  Hertzian  theory  was 
previously  used  by  Crook  [28]  In  a  study  of  impacts  between 
metal  bodies.  The  experimental  results  from  [14]  and 
present  study  show  that  the  value  of  q  varies  from  1.5  to 
2.5.  Local  plastic  deformation,  anisotropic  properties  of 
composite  material  and  unloading  rate  are  all  possible 
causes  for  this  deviation.  Obviously,  an  analytical  study 
to  determine  the  value  of  q  as  function  of  aforementioned 
factors  Is  impracticable.  Since  the  purpose  of  this  study 
is  to  establish  a  contact  law  that  can  be  used  in  the 
analysis  of  impact,  the  validity  of  this  law  must  be 
verified  from  impact  experiment.  This  will  be  investigated 
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in  the  next  chapter. 

From  Equation  (3-3)  or  (3-4),  it  can  be  seen  that  oio 
plays  an  essential  role  in  the  unloading  law  and  hence  the 
value  of  it  must  be  estimated  accurately.  Both  of  Equation 
(3-7)  used  by  Yang  [14]  and  Equation  (3-8)  used  In  this 
study  for  calculating  olq  were  obtained  experimentally,  in 
which  o!cr  arid  (Xp  are  considered  to  be  material  constants  and 
were  determined  using  do  and  am  from  test  data.  However,  it 
was  pointed  out  in  [14]  that  the  values  of  ag  might  not  be 
the  true  permanent  indentations.  They  were  the  values  which 
could  make  the  power  law  given  by  Equation  (3-4)  fit  the 
total  data  under  the  unloading  path.  In  fact,  the  load 
corresponding  to  the  value  of  a^r  -  3.16x10“®  in.  obtained 
in  [14]  is  about  200  lb.  for  0.5  in.  indenter,  which  is 
apparently  too  high.  The  value  of  ttp  =  6.564x10  ^  In. 
obtained  in  this  study,  which  corresponds  to  about  20  1 b  of 
loading,  seems  more  reasonable  as  a  critical  value  in 
indentation.  For  comparison,  the  relations  between 
unloading  rigidity  s  and  maximum  indentation  am  using 
Equation  (3-7)  with  a^r  =  3.16x10"®  In.  and  Equation  (3-8) 
with  ap  =  6.564x10"'*  in.,  respectively,  are  plotted  in 
Figure  3.19.  It  is  interesting  to  see  that  these  two 
equations  give  almost  the  same  values  of  s  up  to  am  =  4x10"® 
in.  which  is  approximately  the  maximum  indentation  before 
failure  could  occur  to  the  specimen.  The  advantage  of  using 
Equation  (3-7)  for  the  formulation  of  the  unloading  law  is 
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that  the  value  of  s  is  constant  for  any  (Xm  once  the  the 
indentation  passes  Wcr.  and  only  one  unloading  test  is 
necessary  to  determine  cXcr  provided  the  load  is  high  enough 
to  produce  permanent  indentations.  The  use  of  Equation  (3- 
9)  needs  performing  many  tests  to  obtain  a  proper  relation 
between  ocq  and  cxm  according  to  Equation  (3-8).  However,  it 
should  be  noted  that  Equation  (3-7)  is  valid  only  if  q  =  5/2 
is  used  In  the  unloading  equation  (3-4),  while  Equation  (3- 
8)  has  no  such  restriction. 
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CHAPTER  4 
IMPACT  EXPERIMENTS 

High  velocity  impacts  usually  result  in  very  small 
contact  time  and  the  material  under  impact  loadings  may 
behave  differently  from  static  contact  due  to  the  strain 
rate  effect.  The  statically  determined  contact  laws 
presented  in  the  previous  chapter  thus  must  be  verified 
experimentally  before  it  can  be  applied  to  the  impact 
analysis.  Wang  [15]  has  conducted  many  impact  experiments 
on  laminated  composite  beams  and  plates  using  spherical 
steel  balls  as  Impacters.  The  strain  response  histories  at 
various  points  on  the  specimens  were  recorded  and  compared 
with  the  finite  element  analysis  with  which  the  contact  laws 
obtained  by  Yang  [14]  was  incorporated.  The  results  showed 
that  the  test  data  agreed  with  the  predictions  using  the 
statical  Indentation  laws  quite  well.  In  this  chapter,  an 
attempt  was  made  to  measure  the  contact  force  directly  so 
that  the  applicability  of  statical  contact  laws  in  impact 
analysis  can  be  further  evaluated. 
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4.1  Experimental  Procedure 

A  6  in.  by  4  in.  laminated  plate  cut  from  a  [0'’/45°/0°/- 
45°/0°]2s  graphite/epoxy  panel  was  used  as  the  impact 
target.  The  0°-direction  was  arranged  to  parallel  the  long 


side  of 

the  plate.  Seven 

stra 

li  n 

gages 

(Micro 

Measurement 

Company 

TYPE  EA-1 3-062 

AQ 

350) 

were 

placed 

at  different 

locations  as  shown  In  Figure 

4.1 

to 

record 

the  dynamic 

strain 

histories.  One 

of 

the 

gages 

was 

placed  on  the 

surface  directly  opposite  to  the  impact  point  to  trigger  the 
oscilloscope.  This  plate  was  hung  with  two  strings  at  two 
corners  to  achieve  the  free  boundary  condition. 

The  projectile  was  made  of  an  impact-force  transducer 
with  a  spherical  steel  cap  of  0.75  inch  in  diameter  glued  on 
the  impact  side  and  a  steel  rod  of  5/8  inch  in  diameter 
glued  on  the  other  side  as  shown  in  Figure  4.2.  It  was  then 
attached  to  a  thin  rod  to  form  a  pendulum  which  could 
produce  impact  velocities  up  to  150in/sec.  The  total  mass 
of  the  projectile  is  0.000181  lb-sec^/ in  . 

The  schematic  diagram  for  this  impact  experimental  set-up 
is  shown  In  Figure  4.3.  Signals  from  gages  and  transducer 
were  amplified  by  a  3A9  Textron ix  amplifier  and  displayed  on 
the  screen  of  an  oscilloscope. 


Figure  4.1  Laminate  dimension  and  strain  gage  locations 


0.75  In  Dia. 


Figure  4.2  Graphical  illustration  of  impact  projectile 
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Figure  4.3  Schemat ical  diagram  for  the  impact 
experimental  set-up 
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4.2  Calibration  of  Impact-Force  Transducer 

The  impact-force  transducer  used  was  Modai  200A05 
marketed  by  PCB  Piezotronics  Inc.  Some  of  It's 
specifications  are  shown  in  Table  4.1  [30].  The  structure 
of  this  transducer  contains  two  thin  quartz  disks  operating 
in  a  thickness  compression  mode  and  sandwiched  between 
hardened  steel  cylindrical  members.  A  built-in  amplifier 
can  reduce  the  high  Impedance  of  the  voltage  from  the  quartz 
element  and  provides  an  output  voltage  which  can  be  read  out 
on  oscilloscope,  recorder,  etc..  The  impact  force  is  then 
computed  using  the  equation, 

F  =  Vp/cp  (4-1 ) 

where  Vp  Is  the  output  voltage  and  Cp  Is  the  sensitivity  of 
the  transducer.  Since  the  value  of  Cp  in  Table  4.1  was 
obtained  under  quasi -static  condition  [30],  it  must  be 
verified  under  Impact  condition  first  so  that  later  the 
results  from  impact  experiment  can  be  correctly  interpreted. 

A  circular  cylindrical  steel  rod  of  2  inch  in  diameter 
and  1.19  Inch  long  hung  on  strings  was  used  as  the  Impact 
target  to  calibrate  the  transducer.  The  acceleration  of  the 
rod  was  measured  by  using  a  Model  302A  accelerometer  which 
was  mounted  on  the  end  of  the  rod  opposite  to  the  impacted 
end  as  shown  in  Figure  4.4.  The  total  weight  of  the  target 
is  1.105  1b. 
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Table  4.1 

Specifications  for  Model  200A05  Impact-Force  Transducer 


Range,  Compression 
(5V  output) 

lb. 

5,000 

Maximum  Compression 

lb. 

10,000 

Resolution  (200  {iM  p-p  noise) 

lb. 

0.2 

St i f f ness 

1  b//ii  n 

100 

Sens i t i V i ty 

mV/ lb 

1  .0 

Resonant  Frequency 
(no  load) 

Hz 

70 , 000 

Rise  Time 

fisec 

10 

Discharge  Time  Constant 
(T.C. ) 

sec 

2,000 

Low-Frequency  (-5%) 

Hz 

0.0003 

L i near i ty , B . F . S . L . 

% 

1 

Output  Impedance 

ohms 

100 

Excitation  (thru  C.C. diode) 

VDC/mA 

+18  to  24/2  to  20 

Temperature  Coefficient 

%/°F 

0.03 

Temperature  Range 

°F 

-100  to  +250 

Shock  (no  load) 

g 

1 0 , 000 
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Using  Equation  (4-1)  and 

a  =  V3/C3  (4-2) 

F  =  ma  (4-3) 

we  obtain 

Cp  =  (4-4) 

where  Vg  and  Ca  are  the  output  voltage  and  the  sensitivity 
of  the  accelerometer,  respectively,  a  is  acceleration  of  the 
target,  and  m  is  the  mass  of  the  target. 

When  impacting  a  metal  projectile  on  a  metal  target  with 
no  pad  on  the  impact  surface,  a  high  frequency  ringing  can 
be  seen  at  the  output  of  the  transducer.  In  order  to  obtain 
smooth  output  curves,  a  soft  pad  was  placed  on  the  Impact 
region  of  the  target  to  eliminate  the  high  frequency 
ringing.  The  cause  of  this  ringing  phenomenon  will  be 
discussed  later.  Typical  output  voltages  of  transducer  and 
accelerometer  read  from  the  oscilloscope  are  shown  in  Figure 
4,5.  Values  of  Vp  were  plotted  vs  the  corresponding  values 
of  Va  taken  from  these  two  curves  at  several  discrete  points 
in  time  and  then  fitted  into  a  straight  line  as  shown  in 
Figure  4.6.  The  slope  of  this  line  represents  the  ratio  of 
Vp/Va  which  is  then  substituted  in  Equation  (4-4)  to 
calculate  the  sensitivity  Cp. 
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Assuming  the  sensitivity  of  the  accelerometer  Ca  is 
correct,  and  using  Equation  (4-4)  and  the  test  data,  the 
average  value  of  Cp  calculated  was  0.494  mV/ lb..  A 
comparison  with  the  value  of  1.0  mV/lb  from  Table  4.1  shows 
that  the  test  result  has  more  than  50%  error.  However, 
since  the  quartz  elements  are  located  at  the  center  of  the 
projectile  while  the  impact  force  is  applied  at  the  end,  we 
were  not  certain  that  the  force  history  picked  up  by  the 
quartz  elements  did  represent  the  real  history  of  the  impact 
force.  The  following  simple  analysis  was  performed  to 
examine  this  uncertainty. 

Consider  a  1  in.  long  steel  rod  with  free-free  boundary 
conditions.  For  a  impulse  loading  given  by 

F(t)  =  Fo  EXP[-(t-T)2/4b2)]  (4-5) 

at  one  end,  the  force  history  at  the  midpoint  of  the  rod, 
Fm(t),  was  computed  and  plotted  in  Figure  4.7  together  with 
the  applied  force  history.  It  should  be  noted  that  the 
values  of  Fq  =  1000  1b.,  r  =  200x10"®  sec.  and  b  =  40x10"® 
sec.  were  chosen  in  Equation  (4-5)  so  that  the  applied 
force  history  is  similar  to  the  experimental  loading 
histroy.  From  Figure  4.7,  it  can  be  seen  that  Fm(t)  is  only 
about  half  of  the  applied  force  F(t).  The  average  ratio  of 
Fm(t)/F(t)  was  obtained  to  be  0.498,  which  is  very  close  to 
the  value  of  Cp  obtained  previously.  The  accelerations  at 
the  two  ends  and  the  midpoint  of  the  rod  were  also 
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calculated  and  plotted  In  Figure  4.8.  It  shows  that  the 
magnitudes  of  acceleration  at  any  position  of  the  rod  have 
virtually  no  difference.  This  indicates  that  the 
accelerometer  did  measure  the  real  acceleration  of  the 
target  while  the  impact-force  transducer  only  picked  up  the 
force  history  at  the  point  of  it’s  own  position.  In  other 
words,  the  wave  motion  in  the  projectile  can  not  be 
neglected,  hence  it  must  be  treated  as  an  elastic  body. 

Repeating  the  previous  analysis  by  changing  the  Impulse 
loading  of  Equation  (4-5)  to 

F(t)  =  FoSin(7rt/b)  (4-6) 

and  letting  Fq  =  1000  lb.  and  b  =  400x10"®  sec.,  we  obtain 
the  force  history  at  the  midpoint  of  the  rod  as  shown  in 
Figure  4.9.  Comparing  Figure  4.9  with  Figure  4.8,  it  is 
clear  that  the  initial  slope  of  the  Impulse  forcing  function 
would  affect  the  amplitude  of  ringing.  The  steeper  the 
Initial  slope  is,  the  higher  the  amplitude  of  ringing  will 
be.  When  impacting  the  steel  projectile  on  graph! te/epoxy 
surface,  this  ringing  phenomenon  was  also  observed. 
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4.3  Finite  Element  Analysis 
4.3.1  Plate  Finite  Element 

A  9-node  Isoparametric  plate  finite  element  (see  Figure 
4.10)  developed  by  Yang  [31]  based  upon  the  laminate  theory 
of  Whitney  and  Pagano  [18]  was  used  to  model  the  dynamic 
motion  of  the  laminated  plate.  At  each  node  there  are  five 
degrees  of  freedom.  Among  them,  u° ,  v°  and  w  are 
displacement  components  of  mid-plane  in  the  x-,y-  and  z- 
dlrectlon,  respectively,  and  <P^  and  <#>y  are  rotations  of  the 
cross-sections  perpendicular  to  the  x—  and  y— axis, 
respectively.  For  symmetric  laminates,  the  flexural 
deformation  Is  uncoupled  from  the  In-plane  extenslonal  and 
shear  deformations,  and  hence,  the  degrees  of  freedom 
corresponding  to  u°  and  v®  can  be  neglected  in  the 
transverse  impact  problem. 

The  Isoparametric  plate  finite  element  Is  developed  using 
the  following  shape  functions: 

For  corner  nodes: 

S.  =  (1/4)(1+§o)(1+77o)(?o+’7o-1  )  +  (1/4)(1-e=')(1-77=)  (4-7) 

For  nodes  at  ?  =  0  and  77  =  ±1  : 

S,  =  (1/2)(1-P)(77o+77=)  ^4-8) 


For  nodes  at  ?  =  ±1  and  77  =  0: 
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The  stiffness  and  mass  matrices  are  obtained  by  numerical 
Integration  using  Gauss  quadrature.  Following  standard 
finite  element  procedures,  the  system  stiffness  matrix  [Kp] 
and  mass  matrix  [Mp]  are  assembled  from  the  element 
matrices.  The  equations  of  motion  are  expressed  in  matrix 
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form  as 

[Mp]{qp}  +  [Kp]{qp}  =  {Pp} 

where 


is  the  force  vector  In  which  F  Is  the  contact  force 
associated  with  the  degree  of  freedom  corresponding  to  the 
w-displacement  at  the  impact  point.  The  subscript  p  in 
Equations  (4-12)  through  (4-15)  denotes  those  are  quantities 
corresponding  to  laminated  plate. 

4.3.2  Modeling  of  Projectile 

In  Section  4.2  we  showed  that  in  order  to  interpret  the 
experimental  transducer  response,  it  is  necessary  to  treat 
the  projectile  as  an  elastic  body.  A  higher  order  rod 
finite  element  developed  by  Yang  and  Sun  [32]  was  used  to 
model  the  projectile.  This  element  has  two  degrees  of 
freedom  at  each  node,  namely  the  axial  displacement  u  and 
it's  first  derivative  9u/9x.  It  has  been  shown  that  this 
higher  order  element  is  far  more  superior  than  the  elements 
with  less  degrees  of  freedom  in  the  analysis  of  dynamic 
problems.  The  displacement  function  Is  taken  as 


u  =  a^  +  a2X  +  asX®  +  a4X® 


(4-16) 
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where  &■,  are  constant  coefficients.  Solving  these 
coefficients  in  terms  of  the  nodal  degrees  of  freedom  and 
substituting  into  Equation  (4-16),  we  obtain 


u  =  {N}Mqp}, 


(4-17) 


where 


{q  }  T  =  {(u),,  Ou/ax)i,  (u)2,  (9u/ax)2} 


(4-18) 


the  vector  of  element  nodal  degrees  of  freedom,  and 


{N}T  =  {fi(x),  f2(x),  f3(x),  f4(x)} 


(4-19) 


in  which 


f, (x)  =  (1  -  x/L)=(1  +  2x/L) 
f2(x)  =  x(1  -  x/L)=^ 
f3(x)  =  xVLMS  -  2x/L) 
f4(x)  =  xVL(x/L  -  1) 

are  shape  functions.  The  subscript  r  in  Equation  (4-17) 
denotes  quantities  corresponding  to  the  rod. 

Using  variational  principle,  the  equations  of  motion  for 
one  element  are  obtaIne(d  as 


[mp]{qp}e  +  [kr]{qr>e  =  {Pr>, 


(4-20) 


where  (Pr^e  the  vector  of  the  generalized  forces 
associated  with  the  nodal  degrees  of  freedom  {qr}e.  [frip]  's 
the  element  mass  matrix  whose  entries  are  given  by 
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=  pAj^f.fjdx  1,J  =  1,2, 3. 4  (4-21) 

and  [krl  Is  the  element  stiffness  matrix  whose  entries  are 
g 1 ven  by  * 


(kr)ij  =  EAj*^f  j  '  f  j  '  dx  1  ,J  =  1,2, 3, 4  (4-22) 

In  Equations  (4-21)  and  (4-22),  p,  E  and  A  are  mass  density. 
Young's  modulus  and  cross-sectional  area  of  the  projectile, 
respectively,  and  L  Is  the  length  of  the  element.  The 
explicit  forms  of  [k^]  and  [m^]  are  given  by 

■  36  3L  -36  3L 

EA  3L  4L2  -3L  -L^ 

[kj  =  -  (4-23) 

30L  -36  -3L  36  -3L 

3L  -L=  -3L  4L2 

and 


■  156 

22L 

54 

-1  3L 

pAL 

22L 

4L2 

13L 

-3L2 

420 

54 

13L 

156 

-22L 

_-13L 

-3L= 

-22L 

r 

to 

(4-24) 


Following  the  usual  manner,  the  system  stiffness  and  mass 
matrices  are  assembled  from  the  element  stiffness  and  mass 
matrices,  and  the  system  equations  of  motion  are  expressed 


as 
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[M,]{qn}  +  [Kj{qr}  =  {Pr} 

where 

{P^}T  =  {F.0,---,0}  (4-26) 

in  which  F  is  the  contact  force  applied  at  the  impacting  end 
of  the  project i 1 e. 


4,4.  Results  and  Discussion 

The  6  in.  by  4  in,  graphite/epoxy  laminate  was  modeled 

by  140  (14  X  10  mesh)  plate  elements  while  the  projectile 
was  modeled  by  20  rod  elements  (see  Figure  4.11).  The  two 
sets  of  equations  (4-14)  and  (4-25)  along  with  the  contact 
laws  given  by  Equations  (3-1),  (3-3)  and  (3-11)  were  solved 
simultaneously.  The  finite  difference  method  with  At  =  0.2 
Msec.  was  used  to  integrate  the  time  variable.  A  coarser 
finite  element  mesh  for  plate  was  used  and  it  was  found  that 
the  present  mesh  yielded  converged  solutions.  A  3- 
Dlmenslonal  analysis  using  112  axlsymmetric  finite  elements 
to  model  the  projectile  was  also  performed,  and  the  results 
showed  the  the  response  at  the  midpoint  of  the  projectile  to 
have  no  significant  difference  comparing  with  the  one 
obtained  by  using  rod  elements. 

An  impact  velocity  of  115  in/sec  was  used  in  the 
experiment.  Figures  4.12-4.17  show  the  strain  response 
histories  at  the  six  locations  picked  up  by  the  strain 


Pt.  of  Impact 


—  3@0.5' 


80  0.375 

-6  in - 


•3@0.5- 


(a)  Plate 


.11  Finite  element  mesh  for  1  ami  anted  plate  and 


102 


gages.  The  results  obtained  using:  the  finite 
methods  and  the  contact  laws  are  also  shown 
figures.  It  Is  evident  that  the  finite  element 
agree  with  the  experimental  data  very  well. 


element 
i n  these 
sol ut I ons 


In  Figure  4.18,  the  experimental  transducer  responses  and 
the  computed  transducer  responses  using  finite  element  are 
plotted  against  time  as  curve  I  and  curve  II,  respectively. 
The  computed  contact  force  history  is  also  plotted  as  curve 
III.  It  can  be  seen  that  the  magnitudes  of  curve  I  and 
curve  II  agree  fairly  well.  The  frequencies  of  ringing  for 
these  two  curves,  however,  are  quite  different.  For  the 
finite  element  results,  the  time  interval  between  two 
consecutive  peaks  of  ringing  Is  approximately  equal  to  the 
time  that  the  longitudinal  stress  wave  needed  to  travel  the 
distance  between  two  ends  of  the  projectile.  This  indicates 
that  the  ringing  is  simply  caused  by  the  transient  wave 
travelling  back  and  forth  in  the  projectile. 


From  Figure  4.18  we  can  see  that  curve  I  has  exact  9 
peaks  in  180  microseconds,  and  the  time  interval  between  two 
consecutive  peaks  is  about  20  microseconds.  It  is  noted 
that  this  transducer  has  a  rise  time  of  10  microseconds  (see 
Table  4.1),  which  is  the  time  it  needs  to  reach  the  maximum 
response.  Any  input  signal  with  period  smaller  than  twice 
of  this  value  will  be  smoothed  out  by  the  transducer,  and 
the  output  signal  may  appear  to  have  lower  frequency.  In 
other  words,  the  period  of  the  output  signal  will  be  at 
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Figure  4.12  Strain  response  history  at  gage  No.1 
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Figure  4.13  Strain  response  history  at  gage  No. 2 
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Figure  4.15  Strain  response  history  at  gage  No. 4 
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Figure  4.17  Strain  response  history  at  gage  No. 6 
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Figure  4.18  Transducer  response  and  contact  force 
histories  from  experimental  and  finite 
element  results 
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least  20  microseconds.  This  might  explain  the  lower 
frequency  of  ringing  in  the  output  voltage  from  the 
transducer . 

The  total  duration  of  contact  for  this  impact  test  is 
about  800  microseconds,  and  multiple  contact  is  also 
observed  from  the  test  data.  Figure  4.19  shows  the 
experimental  transducer  responses  and  the  computed 
transducer  responses  up  to  800  microseconds.  Although  these 
two  results  do  not  matched  very  well  after  the  end  of  the 
first  contact,  it  is  evident  that  the  finite  element 
analysis  does  predict  the  multiple  contact  phenomenon,  and 
the  calculated  total  duration  of  contact  is  also 
approximately  the  same  as  the  test  result. 

Figure  4.20  presents  a  number  of  deformed  configurations 
of  the  laminated  plate  after  impact.  It  is  seen  that  at  the 
point  of  impact,  there  is  a  strong  discontinuity  in  slope  of 
the  transverse  displacement  indicating  the  presence  of  a 
significant  transverse  shear  deformation. 
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Figure  4.19  Transducer  response  histories  from 

experimental  and  finite  element  results  up 
to  800  microseconds 


Figure  4.20  Deformed  configurations  of  laminated  plate 
after  impact 
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CHAPTER  5 

SUMMARY  AND  CONCLUSION 

The  laminate  theory  developed  by  Whitney  and  Pagano  was 
employed  for  studies  of  harmonic  wave  and  propagation  of 
wave  front  in  a  [0V45V0V-45V0°]  2s  graphite/epoxy 
laminate.  The  dispersion  properties  of  flexural  waves  were 
investigated.  The  wave  front  surface  was  constructed  using 
ray  theory.  It  was  shown  that  due  to  the  anisotropic 
properties  of  composite  laminate,  the  transient  wave  would 
propagate  with  different  velocities  in  different  directions. 
The  growth  and  decay  of  the  wave  front  strength  were  also 
discussed. 

The  contact  laws  between  0.5  inch  and  0.75  inch  spherical 
steel  Indenters  and  the  graph! te/epoxy  laminate  were 
determined  experimentally  by  means  of  a  statical  indentation 
test.  Loading,  unloading  and  reloading  curves  were  fitted 
into  power  equations.  Linear  relation  was  found  between  the 
permanent  indentation  and  the  maximum  indentation  at 
unloading,  which  is  seen  to  be  independent  of  the  size  of 
indenters.  This  relation  was  then  used  to  determine  the 
coefficient  of  the  unloading  law.  It  was  demonstrated  that 
there  was  no  need  to  perform  reloading  experiments  once  the 
loading  and  unloading  laws  were  established. 


Test  results 
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showed  loading  and  reloading  curves  followed  the  power  laws 
with  power  indices  of  1.5  very  well,  while  the  power  indices 
for  unloading  curves  varied  from  1.5  to  2.5. 

The  statically  determined  contact  i aws  were  incorporated 
into  an  existing  9-node  isoparametric  plate  finite  eiement 
program  to  study  the  dynamic  response  of  a  graph! te/epoxy 
laminated  plate  subjected  to  impact  of  a  hard  object.  An 
impact  experiment  was  conducted  to  verify  the  validity  of 
statical  contact  laws  in  the  dynamical  impact  analysis.  It 
was  shown  that  the  strain  responses  predicted  using  the 
finite  element  method  agreed  with  the  test  results  very 
well.  The  contact  force  history  of  the  impact  test  was 
measured  by  an  impact-force  transducer,  which  was  also  seen 
to  match  the  finite  element  result  In  magnitude  as  well  as 

contact  duration. 

The  indentation  tests  have  been  used  ever  since  the 
beginning  of  the  century  to  determine  the  static  and  dynamic 
hardnesses  of  metals  in  terms  of  the  applied  loading,  the 
size  of  the  indenter,  and  the  chordal  diameter  of  the 
permanent  indentation  [33].  If  similar  systematic 
indentation  tests  are  performed  on  the  laminated  composite 
materials,  then  the  relations  between  contact  coefficients 
and  the  sizes  of  the  indenters  could  be  determined  more 
rigorously,  and  the  usefulness  of  the  contact  laws  could  be 


further  extended. 
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As  the  verification  of  the  contact  laws  has  been  limited 
to  low  velocity  Impacts  In  this  study,  their  accuracy  under 
high  velocity  Impact  conditions  Is  not  clear.  Besides  the 
contact  behavior  which  may  be  significantly  different  from 
the  static  one,  the  damage  Induced  by  waves  could  be  quite 
extensive  which  needs  to  be  Included  In  the  analysis.  While 
the  present  study  tried  to  establish  experimentally  contact 
laws  which  can  be  used  In  the  analysis  of  low  velocity 
Impact,  the  damage  of  laminate  due  to  Impact  loading  has  not 
been  discussed.  It  is  apparent  that  more  work  needs  to  be 
done  so  that  the  failure  mechanism  In  laminated  composites 
due  to  Impact  can  be  better  understood.  Stress  waves 
propagating  in  thickness  direction,  which  may  be  responsible 
for  the  del  ami  nation  of  laminates.  Is  one  of  the  Important 
subjects  that  should  be  Investigated.  Strength  and  fatigue 
life  degradations  of  laminates  after  Impact,  which  have  been 
examined  briefly  by  Wang  [15],  also  need  more  extensive 
study . 
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APPENDIX 

COMPUTER  PROGRAM  AND  USER  INSTRUCTIONS 

The  computer  program  used  in  this  research  was  written 
following  the  program  by  Professor  R.  L.  Taylor  [34]  with 
some  necessary  modification  in  order  to  solve  the  impact 
problems  of  laminated  plates.  A  brief  Instruction  of  the 
input  data  for  solving  the  impact  problem  specified  in 
Chapter  4  of  this  report  is  given  in  this  apppendix.  The 
detailed  descriptions  of  data  input  as  well  as  the  macro 
instructions  for  solving  various  types  of  problems  can  be 
found  in  [34].  The  listing  of  input  is  shown  at  the  end  of 
this  appendix,  followed  by  the  listing  of  program. 

I.  Title  and  control  information: 

1.  Title  card-Format(20A4) 

Col umns  Descr i pt i on 

1-4  Must  contain  FECM 

5- 80  Alphanumeric  information  to  be  printed  with 

output  as  page  header. 

2.  Control  information  card-Format (615) 

Col umns  Descr i pt I  on 

1-5  Number  of  nodes  (NUMNP) 

6- 10  Number  of  elements  (NUMEL) 
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11-15  Number  of  layers  (LAYER) 

16-20  Spatial  dimension  (NDM) 

21-25  Number  of  unknowns  per  node  (NDF) 

26-30  Number  of  nodes  per  element  (NEN) 

II.  Mesh  and  initial  information: 

The  input  of  each  segment  in  this  part  of  data  is 
controlled  by  the  alphanumeric  value  of  macros,  which  must 
be  followed  immediately  by  the  appropriate  data.  Except  for 
the  END  card  which  must  be  the  last  card  of  this  part,  the 
data  segemnts  can  be  in  any  order.  Each  segment  is 
terminated  with  blank  card(s).  The  meaning  of  each  macro  is 
given  by  the  following: 

Macro  Data  to  be  input 

COOR  Coordinate  data 

ELEM  Element  data 

BOUN  Boundary  condition  data 

MATE  Material  data 

ROD  Initial  condition  of  the  project i 1 e 

EXPE  Experimental  indentation  laws  data 
END  Must  be  the  last  card  of  this  part,  terminates 

mesh  and  initial  information  input. 

1.  Coordinate  data-Format(2I5,2F10.0) 

Col umns  Descr i pt i on 

1-5  Nodal  number 

6-10  Generation  increment 


11-20  X-coordinate 
21-30  Y-coordinate 

2.  Element  data-Format(1 1 15) 

Col umns  Descr i pt I  on 

1-5  Element  number 
6-10  Node  1  number 
11-15  Node  2  number 
etc . 

46-50  Node  9  number 
51-55  Generation  increment 

3.  Boundary  condition  data-Format (715) 

Col umns  Descr i pt i on 

1 -5  Node  number 
6-10  Generation  increment 
11-15  DOF  1  boundary  code 

1 6-20  DOF  2  boundary  code 

21-25  DOF  3  boundary  code 

26-30  DOF  4  boundary  code 

31-35  DOF  5  boundary  code 

4.  Initial  condition  of  the  project! le-Format(2I5,F10.0) 
Col umns  Descr i pt i on 

1-5  The  node  at  which  the  projectile  hits 
6-10  DOF  corresponding  to  the  direction  of  impact 
11-20  Initial  impact  velocity 
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5.  Experimental  Indentation  laws  data-Format(4F10.0) 

Co I umns  Descr I pt I  on 

1-10  Contact  coefficient  k 
11-20  Critical  Indentation  Wp 

21-30  Constant  Sp  of  Equation  3-9 

31-40  Power  index  q  of  the  unloading  law 

6.  Material  data 

Card  1-format(3!5,F10.0) 

Col umns  Descr I pt I  on 

1-5  Order  of  Gauss  quadrature  for  the  numerical 
integration  of  the  bending  energy 
6-10  Order  of  Gauss  quadrature  for  the  numerical 
integration  of  the  transverse  shear  energy 
11-15  Order  of  Gauss  quadrature  for  strain  outputs 

at  Gauss  points  If  >0 

at  nodal  points  if  <0 

16-25  Total  thickness  of  the  laminate 

Card  2-Format (7F1 0 . 0) 

Columns  Description 


1-10 

Mass  density 

11-20 

Poisson's  ratio  Vi2 

21-30 

Longitudinal  Young's  modulus  E 

31-40 

Transverse  Young's  modulus 

41-50 

Shear  modulus  0^2 

11-20 

Shear  modulus  G^ 3 

11-20 

Shear  modulus  G23 
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Card  3,4: 


Format ( I5,F5.0,F10.0) 


Col umns  Descr i pt I  on 

1-5  Layer  number 

6-10  Fiber  angle 

11-20  Thickness  of  the  layer 


III.  Macro  instructions: 


The  first  instruction  must  be  a  card  with  MACR  in  columns 

1  to  4.  The  macro  instructions  needed  to  solve  the  problem 

specified  in  Chepter  4  of  this  report  are  shown  in  the 

listing  of  input.  Cards  must  be  input  in  the  precise  order. 

The  following  is  the  explanation  of  each  macro: 

Columns  Columns  Columns 

1-4  5-10  11-15  Descr i pt i on 


LMAS 


LOOP 


TIME 


RODP 


DISP 


Lumped  mass  formulation 
Set  time  increment  to  value  V 
Execute  N  times  the  instructions 
between  this  macro  and  macro  NEXT 
Advance  time  by  DT  value 
Integration  of  the  equations  of 
motion  using  the  finite  difference 
method.  Contact  force,  indentation 
and  element  strain  will  be  stored 
stored  every  N  steps  in  loop 
Nodal  displacements  will  be  stored 
every  N  steps  in  loop 
End  of  loop  instructions 


NEXT 
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END  End  of  macro  program  Instructions 

IV.  Termination  of  program  execution 

A  card  with  STOP  in  columns  1  to  4  must  be  supplied  at 
the  end  of  the  input  data  in  order  to  properly  terminate  the 
execut i on . 

The  values  of  contact  force,  indentation,  element  strain, 
nodal  displacement  and  the  response  of  the  projectile  at 
each  requested  output  time  step  are  stored  In  program  files 
which  can  be  saved  (say,  copy  to  a  magnetic  tape)  at  the  end 
of  execution.  Three  program  files,  i.e.;  tapeS,  tapeS  and 
tape9  ar©  used  for  data  saving: 

TapeS:  Nodal  displacement  -  Format (6E1 2.4) 

Nodal  displacements,  from  node  1  to  node  NUMNP,  are  saved 
on  tapes  at  each  requested  output  time  step  according  to  the 
format . 

TapeS:  Element  strain  -  Format (21 6, 5E1 2.4) 

Element  strains,  from  element  1  to  element  NUMEL,  and 
then  from  node  1  to  node  NEN  of  each  element,  are  saved  on 
on  tapes  at  each  requested  output  time  step. 

Co  1 umns  Data  saved 

1-6  Element  number 

7-12  Node  number  of  element 

1  S-24  Bending  strain 


25-S6  Bending  strain  /Cy 
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37-48  Bending  strain 

49-60  Transverse  shearing  strain 

49-60  Transverse  shearing  strain 

Tape9:  Contact  force,  indentation  and  the  response  of  the 
projectile  -  Format (6E1 2 . 4) 

The  following  information  is  saved  on  tape9  at  each 
requested  output  time  step: 

Columns  Data  saved 


1-12 

Contact  force 

13-24 

I ndentat i on 

25-36 

'Transducer ' 

response  (see  Chapter 

4) 

37-48 

Displacement 

of 

the  projectile  at 

the 

impacted 

end 

37-48 

Velocity  of 

the 

projectile  at  the 

impacted  end 

37-48 

Accel erat i on 

of 

the  projectile  at 

the 

impacted 

end 
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LISTING  OF  INPUT  DATA 


FECM  *»LOW  UELOCITY  IMPACT  OF  LAMINATED  PLATE»« 
G09  140  20  2  5  9 


1 

1 

0.0 

0.0000 

7 

1 

1.5 

0.0000 

23 

1 

4.5 

0.0000 

29 

0 

G.O 

0.0000 

30 

1 

0.0 

0.2500 

3B 

1 

1.5 

0.2500 

52 

1 

4.5 

0.2500 

58 

0 

G.O 

0.2500 

59 

1 

0.0 

0.5000 

G5 

1 

1.5 

0.5000 

81 

1 

4.5 

0.5000 

87 

0 

G.O 

0.5000 

88 

1 

0.0 

0.6875 

94 

1 

1.5 

0.6875 

no 

1 

4.5 

0.B875 

IIG 

0 

6.0 

0.G875 

117 

1 

0.0 

0.8750 

123 

1 

1.5 

0.8750 

139 

1 

4.5 

0.8750 

145 

0 

G.O 

0.8750 

14G 

1 

0.0 

1.0B25 

152 

1 

1.5 

1.0G25 

188 

1 

4.5 

1 . 0G25 

174 

0 

6.0 

1.0G25 

175 

1 

0.0 

1.2500 

181 

1 

1.5 

1.2500 

197 

1 

4.5 

1.2500 

203 

0 

G.O 

1.2500 

204 

1 

OoO 

1.4375 

210 

1 

1.5 

1.4375 

22B 

1 

4.5 

1.4375 

232 

0 

G.O 

1.4375 

233 

1 

0.0 

1.G250 

239 

1 

1.5 

1.G250 

255 

1 

4.5 

1.G250 

2B1 

0 

G.O 

1.G250 

2G2 

1 

0.0 

1.8125 

2G8 

1 

1.5 

1.8125 

284 

1 

4.5 

1.8125 

290 

0 

G.O 

1.8125 

291 

1 

0.0 

2.0000 

297 

1 

1.5 

2.0000 

313 

1 

4.5 

2.0000 

319 

0 

G.O 

2.0000 

320 

1 

0.0 

2.1875 

32G 

1 

1.5 

2.1875 

342 

1 

4.5 

2.1875 

348 

0 

G.O 

2.1875 

349 

1 

0.0 

2.3750 

355 

1 

1.5 

2.3750 

371 

1 

4.5 

2.3750 

377 

0 

G.O 

2.3750 

378 

1 

0.0 

2.5G25 

384 

1 

1.5 

2.5B25 

400 

1 

4.5 

2.5G25 

40B 

0 

G.O 

2.5625 

407 

1 

0.0 

2.7500 

413 

1 

1.5 

2.7500 

423 

1 

4.5 

2.7500 

435 

0 

G.O 

2.7500 

43G 

1 

0,0 

2.9375 

442 

1 

1.5 

2.9375 

127 


458  1 

4.5 

4B4  0 

6.0 

4G5  1 

0.0 

471  1 

1.5 

487  1 

4.5 

433  0 

6.0 

494  1 

0.0 

500  1 

1.5 

51G  1 

4.5 

522  0 

6.0 

523  1 

0.0 

529  1 

1.5 

545  1 

4.5 

551  0 

6.0 

552  1 

0.0 

558  1 

1.5 

574  1 

4.5 

580  0 

6.0 

581  1 

0.0 

587  1 

1.5 

B03  1 

4.5 

B09  0 

6.0 

ELEN 

1  1 

3  G1 

15  59 

B1  119 

23  117 

119  177 

43  175 

177  235 

57  233 

235  293 

71  291 

293  351 

85  343 

351  409 

93  407 

409  457 

113  465 

4B7  525 

127  523 

525  583 

BOUN 

1  1 

-1  -1 

B09  0 

1  1 

ROD 

505  3 

115.0 

EXPE 

1912000. 

0.000B5G4 

NATE 

3  3 

-3 

0.000148 

0.3 

1  0. 

0.0053 

2  45. 

0.0053 

3  0. 

0.0053 

4  -45. 

0.0053 

5  0. 

0.0053 

G  0. 

0.0053 

7  45. 

0.0053 

8  0. 

0.0053 

9  “45. 

0.0053 

10  0. 

0.0053 

11  0. 

0.0053 

12  -45. 

0.0053 

13  Oo 

0.0053 

14  45. 

0.0053 

15  0. 

0.0053 

16  0. 

0.0053 

17  -45. 

0.0053 

18  0. 

0.0053 

13  45. 

0.0053 

20  0. 

0.0053 

END 


2.9375 

2.9375 

3.1250 

3.1250 

3.1250 

3.1250 

3.3125 

3.3125 

3.3125 

3.3125 

3.5000 

3.5000 

3.5000 

3.5000 

3.7500 

3.7500 

3.7500 

3.7500 

4.0000 

4.0000 

4.0000 

4.0000 


59 

2 

32 

GO 

117 

GO 

90 

118 

175 

118 

148 

176 

233 

176 

206 

234 

291 

234 

264 

292 

349 

292 

322 

350 

407 

350 

380 

408 

465 

408 

438 

466 

523 

466 

496 

524 

581 

524 

554 

582 

0 

0 

0 

0 

0 

0 

0.094  2.0 

.lOG 

17500000.  1150000. 


30  31 

88  89 

146  147 

204  205 

262  263 

320  321 

378  379 

436  437 

494  495 

552  553 


800000. 


2 

2 

2 

2 

2 

2 

2 

2 

2 

2 


800000.  800000. 
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MACR 

LMAS 

DT 

LOOP 

TIME 

RODP 

DISP 

NEXT 

EMD 

STOP 


.2E-G 

10 


5 

5 
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LISTING  OF  PROGRAM 


PROGRAM  MAIM ( INPUT » OUTPUT  *  TAPE5=INPUT » T APEG-OUTPUT>  TAPES*  TAPES* 
1  TAPES* TAPES) 

MAIN  PROGRAM 
LOGICAL  PCOMP 

COMMON  ^TDATA/  2?HEAD(S0 ),NUMNP*NUMEL* LAYER* NEQ.IPR 
COMMON  /LABELS/  PDIS(B)*A(B)*BC(3)*DI(G)*CD(3)*FD(3) 

COMMON  /LODATA/  NDF*NDM*NEN*NST*NKM 
COMMON  /PARATS/  NPARC 14) . NEND 
DIMENSION  TITL(aO)*WD(3) 

COMMON  G(3S000) 

DIMENSION  MOSOOO) 

EQUIUALENCE  (G(1)*N(1)) 

MAX=33000 
WD(1)=4HFECM 
UD(a)=4HMACR 
WD(3)=4HST0P 
SSS  READ(5.1000)  TITL 

IF(PCOMP(TITL(l)*MD(l)))  GO  TO  100 
IF(PCOMP(TITL(l)*WD(a)))  GO  TO  SOO 
IF(PC0MP(TITL(1).WD(3)))  STOP 
GO  TO  339 

100  DO  101  1=1,30 

101  HEAD(I)=TITL(I) 

READ(5, 1001)  NUMNP*NUMEL* LAYER* NDM*NDF* MEN 
WRITE(6*3000)  HEAD*NUMNP*  NUMEL  *  LAYER  *  NDM  *  NDF  *  NEN 
PDIS(3)=A(NDM) 

NST=NEN«-NDF 
DO  110  1=1*14 
no  NPAR(I)=1 
NPARC 1)=1 

NP  AR  ( a )  =NPAR  ( 1 )  +3-»NST»IPR 
NP  AR  ( 3 )  =NPAR  ( a ) +NDM^fNEN»  I  PR 
NP AR ( 4 ) =NPAR ( 3 ) +NST 
NP AR ( 5 ) =NPAR ( 4 ) +NST « I PR 
NPAR  ( G )  =NPAR  ( 5 )  +NEN-:f  NUMEL 
NPARC?)  =NPAR  C  G )  +NDF«-NUMNP 
NPAR  C  8 ) =NPAR  C  7 ) +NDM«NUMNP» IPR 
NP AR  C  9 ) =NPAR  C  8 ) +NDF »NUMNP»I PR 
NP AR  C 1 0 ) =NP AR  C  9 ) +NDF*NUMNP 
CALL  SETMEMCNPARC9)) 

CALL  PZER0CGC1),NPARC9))  ^ 

CALL  PMESH  C  M  C  NP  AR  C  3 ) )  *  G  C  NPAR  C  a  ) )  *  M  C  NPAR  C  5  ) )  *  ^  (  MPAR  (G  )  )  * 

1  G  C  NPAR  C  7 ) )  *  G  C  NPAR  C  8 ) )  *  M  C  NP  AR  C  9 ) )  *  NDF  *  NDM  *  NEN  *  NKM ) 

NPARC 10)=NPARC3)+NEQ 

NPARC 11 )=NPARC 10 )+NDF«NUMNP»IPR 

NEND=NPARC 1 1 )+NEQ*IPR 

NE=NEND 

CALL  SETMEMCNE) 

CALL  PZERO  C  G  C  NP AR  C 1 0 ) ) *  NE-NPAR  CIO)) 

GO  TO  S93 

aoo  CALL  PMACRCGCNPARCl))*GCNPARCa)),MCNPARC3)),G(NPARC4)), 

1  MCNPARC5)),MCNPARCG)),GCNPARC7))*GCNPARC8))*MCNPARC9)), 

a  G  C  NPAR  C 1 0 ) ) *  G  C  NPAR  C 1 1 ) ) *  G  C  NE ) *  NDF  *  NDM*  NEN*  NST ) 

CALL  PZERO CG* MAX) 
rn  Tfi 

1000  FORi1AT(20A4) 

1001  F0RMm'(lGI5) 

2000  FORMOT(1H1.20A4//  ^  „  k.  . 


5X9?^C  0  N  T  R  0  L 


0  R  M  A  T  I  0  N  S?^// 


10X.35HNUMBER  OF  NODAL  POINTS 
lOX^SSHNLiMBER  OF  ELEMENTS 
10Xp35HNIJMBER  OF  MATERIAL  LAVERS 
10X9  35HDII1EN5I0N  OF  COORDINATE  SPACE 


=,16/ 

=,16/ 

=,16/ 

=,16/ 


iOX,35HDEGREES  OF  FREEDOM  FOR  EACH  NODE  =,16/ 
10X,35HN0DES  PER  ELEMENT  (MAXIMUM)  =, IS) 


MAIN  1 
MAIN  2 
MAIN  3 
MAIN  4 
MAIN  5 
MAIN  6 
MAIN  7 
MAIN  8 
MAIN  9 
MAIN  10 
MAIN  11 
MAIN  12 
MAIN  13 
MAIN  14 
MAIN  15 
MAIN  16 
MAIN  17 
MAIN  18 
MAIN  19 
MAIN  20 
MAIN  21 
MAIN  22 
MAIN  23 
MAIN  24 
MAIN  25 
MAIN  26 
MAIN  27 
MAIN  28 
MAIN  29 
MAIN  30 
MAIN  31 
MAIN  32 
MAIN  33 
MAIN  34 
MAIN  35 
MAIN  36 
MAIN  37 
MAIN  38 
MAIN  39 
MAIN  40 
MAIN  41 
MAIN  42 
MAIN  43 
MAIN  44 
MAIN  45 
MAIN  46 
MAIN  47 
MAIN  48 
MAIN  49 
MAIN  50 
MAIN  51 
MAIN  52 
MAIN  53 
MAIN  54 
MAIN  55 
MAIN  56 
MAIN  57 
MAIN  58 
MAIN  59 
MAIN  60 
MAIN  61 
MAIN  62 
MAIN  63 
MAIN  64 
MAIN  65 
MAIN  66 
MAIN  67 
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C  •  •  •  • 


BLOCK  DATA 
BLOCK  DATA 

COMMON  /CTDATA/  O.HEAD(EO)»NLJMMP.NUMEL»LAYER»NEQ. IPR 
COMMON  /LABELS/  PDIS(6)9A(B)»BC(2)9DI(B)?CD(3)»FD(3) 

DATA  O/lHl/.IPR/l/ 

DATA  PDIS/4H(I10»EH9  94HF13.»4H4»  94HGE13.4H.4)  / 

DATA  A/EHj  1 9  EH>  Es.  EHp  3?  eh?  4?  EH?  5?  EH?  G/ 

DATA  BC/4H  B.C?EH.  / 

DATA  DI/4H  DIS?EHPL?4H  UEL?EH0C?4H  ACC?EHEL/ 

DATA  CD/4H  COO? 4HRDIN? 4HATES/ 

DATA  FD/4H  FOR? 4HCE/D? 4HISPL/ 

END 

SUBROUTINE  PMftCR(UL?XL?LD?P? IX? ID?X?F? JDIAG?DR?B?CT?NDF?NDM? 
NEN?NST) 

MACRO  INSTRUCTION  ROUTINE 
LOGICAL  PCOMP 
COMMON  G(l) 

DIMENSION  M(l) 

EQUIUALENCE  (G(1)?M(1))  „ 

COMMON  /CTDATA/  0? HEAD(EO) ? NUMNP? NUMEL? LAYER? NEQ? IPR 
COMMON  /PROLOD/  PROP 

COMMON  /TMDATA/  TIME? DT? DDT? FORCE? ALPHA 
COMMON  /ISNIDX/  ISU 
COMMON  /PARATS/  NPAR(14)?NEND 
COMMON  /RODATA/  UR?IQ?NDS 

DIMENSION  UL(1)?XL(1)?LD(1)?P(1)?IX(1)? ID(1)?X(1)?F(1)? 

^  JDIAG(1)?DR(1)?B(1) 

DIMENSION  WD(b)?CT(4?  lG)?LUEO) 

DATA  WD/4HL00P?4HNEXT?4HDT  ,4HPR0P?4HLMAS?4HR0DP? 

1  4HSTRE?4HDISP?4HCHEC/ 

DATA  NWD/S/.ENDM/4HEND  / 

INITIALIZATION 
DT  =0,0 
PROP  =1.0 
TIME  =  0.0 
NNEQ  =  NDF^fNUMNP 
NPLD  =  0 
FORCE=  0, 

ALPHA=  0. 

WRITE(G?E001)  0?HEAD 
LL  =  1 
LMAX  =  IG 

CALL  SETMEM(NEND+LMAX«4»IPR) 

CT(1?1)  =  ND(1) 

CT(3?1)  =  1.0 
I  LL  =  LL  +  1 
IFCLL.LT.LMAX)  GO  TO  110 
LMAX  =  LMAX  +  IG 
CALL  SETMEM(NEND+LMAX-»4’'aPR) 

I  READ(5?1000)  (CT(J?LL)? J=l?4) 

WRITE(G?2000)  (CT( J?LL)? J=l?4) 

IF(.NOT.PCOMP(CT(l?LL)?ENDM))  GO  TO  100 
CT(1?LL)  =  WD(E) 

NEND  =  MEND  +LMAX-»4JsIPR 
LX  =  LL  -  1 
DO  230  L=1?LX 

IF(.NOT.PCOMP(CT(l?L)?WD(l)))  GO  TO  E30 
J  =  1 
K  =  L  +  1 
DO  ElO  I=1<?LL 

IF(PCOMP(CT(l?I)?WD(l)))  J  =  J  +  1 
IF(J  .GT,  9)  GO  TO  401 
IF(PC0MP(CT(l9l)?WD(S)))  J  =  J  -  1 
I  IF(J.EQ.O)  GO  TO  220 
GO  TO  400 
I  CT(4.I)  =  L 
CT(4?L)  =  I 
I  CONTINUE 
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j  =  0 

DO  240  L=1,LL 

IF(PCOMP(CT(l»L)»WD(l)))  J  =  J  +  1 
240  IF(PC0MP(CT(1»L),I4D(2)))  J  =  J  -  1 
IF(J-NE.O)  GO  TO  400 
LU  =  0 
L  =  1 

299  DO  300  J=l,mD 

300  IF(PCOriP(CT(lpL)jWD(J)))  GO  TO  310 
GO  TO  330 

310  I  =  L  -  1 

GO  TO  (1p29394.5pG,7,8,9), J 
C _  BET  LOOP  START  INDICATORS 

1  LU  =  LU  +  1 
LX  =  CT(4,L) 

LUE(LU)  =  LX 
CTCS^LX)  =  Ip 
GO  TO  330 

C....  LOOP  TERMINATOR  CONTROL 

2  N  =  CT(4.L) 

CT(3.L)  =  CT(3.L)  +  1.0 
IF(CT(3*L)pGToCT(3.N))  LU  =  LU  -  1 
IF(CT(3,L)pLE.CT(39N))  L  =  N 
GO  TO  330 

C....  SET  TIME  INCREMENT 

3  DT  =  CT(39L) 


DDT=  DT-DT 
GO  TO  330 

C...P  INPUT  PROPORTIONAL  LOAD  TABLE 

4  NPLD  =  CT(3,L) 

PROP  =  PR0PLD(0p,NPLD) 

GO  TO  330 

C _  FORM  LUMPED  MASS  MATRIX 

5  ISW=3 
CALL  KMLIB 
GO  TO  330 

Co...  IMPACT 
G  NDS=CT(3,L) 

IF(NDSpEQoO)  NDS=1 
CALL  RODIPCT 
GO  TO  330 

C....  PRINT  STRESS/STRAIN  UALUE 
7  ISW=4 


LX  =  LUECLU) 

IF(AM0D(CT(39LX).AMAXl(CT(3,L)a.)))  330,71,330 
71  CALL  FSTREA(UL,XL,LD,P, IX, IDpX^F, JDIAG,DR,B,NDF,NDM,NEN,NST,NNEQ) 
GO  TO  330 

C _  PRINT  DISPLACEMENTS 

8  LX  =  LUECLU) 

irCAMODCCTOj-LX),  AMAX1(CT(3j.L),  1. )))  330,81,330 
81  CALL  FRTDISCUL, ID,X9B,F,DR,NDM,NDF) 


GO  TO  330 

Cpopp  check 

9  WRITE(G,5001)  NEND, JDIAGCNEQ) 


330 


RETURN 

L=L-:-l 

IFCLpGTpLL)  return 


GO  TO  233 

Ccpp.  PRINT  ERROR  FORMATS 
400  RRITE(G,4000) 


RETURN 

401  WRITE(6,4001) 

RETURN 

Cpoop  INPUT/OUTPUT  FORMATS 
1000  FORMAT(A4,1X,A4,1X,2F5pO) 

2000  FOF^MATCiOX,  A4,  IX,  A4, 1X,2G15.5) 

2001  FORMAT(A1,20A4//,5X, ISHMACRO  INSTRUCTI0NS//5X, 15HMACR0  STATEMENT 
-  ,5X,iOMUARIABLE  1 , 5X, lOHUARIABLE  2) 

4000  FORMAT (SX,4GH--^v*PMACR  ERROR  01--*  UNBALANCED  LOOP/nEXT  MACROS  ) 

4001  FORMAT (5Xp45H-«-PMACR  ERROE  02->^  LOOPS  NESTED  DEEPER  THAN  8) 
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PMAC  73 
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PMAC  84 
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PMAC  86 
PMAC  87 
PMAC  88 
PMAC  89 
PMAC  90 
PMAC  91 
PMAC  92 
PMAC  93 
PMAC  94 
PMAC  95 
PMAC  SB 
PMAC  97 
PMAC  98 
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PMAC112 
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PMAC123 
PMAC124 
PMAC125 
PMAC126 


5001  F0RMAT(1H1p////5Xp32HCHECK  MESH  DATA  AND  MEMORY  SPACE// 

10X,12H  MEND  =,  I10//10X,  12HJDIAG(NEQ)  =,110) 

END 

C 

SUBROUTINE  PZERO(U,NN) 

ZERO  REAL  ARRAY 
DIMENSION  U(NN) 

DO  100  N=1,NN 
100  U(N)  =  0„0 
RETURN 
END 
C 

SUBROUTINE  SETMEMCJ) 

MONITOR  AUAIABLE  MEMORY  IN  BLANK  COMMON 
COMMON  /PRSIZE/  MAX 
K  =  J 

IF(K,LE.MAX)  RETURN 
WRITECG, 1000)  K,MAX 
STOP 

1000  FORMAT (5X 9  49hLi'"'^SETMEM  ERROR  01^^  INSUFFICIENT  STORAGE  IN  BLANK, 
-  8H  COMMON  //ITXpIIHREQUIRED  =, I8/17Xp IIHAUAILABLE  =,18) 

END 

C 

LOGICAL  FUNCTION  PCOMP(ApB) 

LOGICAL  COMPARISON 
IFCA-B)  IOpOOpIO 
10  PCGMP  =  -  FALSE o 
RETURN 

20  PCOMP  =  cTRUEo 
RETURN 
END 
C 

SUBROUTINE  ACTCOLC A, B, JDIAG, NEQ, AFAC, BACK, ISS) 

ACTIUE  COLUMN  PROFILE  SYMMETRIC  EQUATION  SOLUER 
LOGICAL  AFAC, BACK, FLAG 
DIMENSION  A(1),B(1), JDIAG(l) 

C..«.  FACTOR  A  TO  UT^D^U,  REDUCE  B 
FLAG=. FALSE. 

JR  =  0 

DO  BOO  J=1pNEQ 
JD  =  JDIAG(J) 

JH  =  JD  -  JR 
IS  =  J  -  JH  +  2 
IF(JH-2)  600,300,100 
100  IFC.NOT.AFAC)  GO  TO  500 
IE  =  J  “  1 
K  =  JR  +  2 
ID  =  JDIAG(IS“i) 

C....  REDUCE  ALL  EQUATIONS  EXCEPT  DIAGONAL 

DO  200  I=IS,IE 
IR  =  ID 
ID  =  JDIAG(I) 

IH  =  MINOCID-IR-IpI-IS+I) 

IF(IH.GT.O)  A(K)=A(K)--D0T(A(K“IH),A(ID-"IH),IH) 

200  K  =  K  +  1 

C.O..  REDUCE  DIGONAL  TERM 

300  IFC.NOT.AFAC)  GO  TO  500 
IR  =  JR  +  1 

IE  =  JD  -  1 
K  =  J  -  JD 
DO  400  I=IR,IE 
ID  =  JDIAG(I<+I) 

IFCACID))  301,400,301 

301  D  =  ACI) 

ACI)  =  A(I)/A(ID) 

ACJD)  =  A(JD)  -  D-^A(I) 

400  CONTINUE 

IF(A(JD))450,450,500 
450  IFCISSoNE.O)  GO  TO  500 
IF (FLAG)  GO  TO  465 
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WRITE(6,460) 

460  FORMAT (//50H»«’ACTCOL  ERROR  01*^^  STIFFNESS  MATRIX  NOT  POSITIUE 
1  8HDEFINITE) 

FLAG=.TRUE. 

465  WRITE(G.46G)  J^ACJD) 

466  F0RMAT(32H  N0NP03ITIUE  PIUOT  FOR  FOLIATION  , I4» 5Xf 7HP0UIT  =» 

-  E20»10) 

C.*..  REDUCE  RHS 

500  IF(BACK)  B(J)  =  B(J)  ~  DOT(A( JR+DpBCIS-I). JH-1) 

GOO  JR  =  JD 

IF (FLAG)  STOP 
IF (.NOT. BACK)  RETURN 
C....  DIUIDED  BV  DIAGONAL  PIUOTS 
DO  700  1=1, NEQ 
ID  =  JDIAG(I) 

IF(A(ID))  650,700,650 
650  B(I)  =  B(I)/A(ID) 

700  CONTINUE 

C....  BACK  SUBSTITUTE 
J  =  NEQ 
JD  =  JDj:AG(J) 

800  D  =  B(J) 

J  =  J  -  1 
IF(J.LE.O)  RETURN 
JR  =  JDIAG(J) 

IF(JD-JR.LE.l)  GO  TO  1000 
IS  =  J  -  JD  +  JR  +  2 
K  =  JR  -  IS  +  1 
DO  300  I=IS,J 

900  B(I)  =  B(I)  -  A(I+K)*D 
1000  JD  =  JR 
GO  TO  800 
END 
C 

SUBROUTINE  ADDSTF(A,S,P, JDIAG,LD,NST,NEL,FLG) 
assemble  global  ARRAYS 
LOGICAL  FLG 

dimension  a ( 1 ) , S ( NST, 1 ) , P ( 1 ) , JDI AG( 1 ) , LD( 1 ) 

DO  200  J=1,NEL 
K  =  LD(J) 

IF(KoEQ.O)  GO  TO  200 
IF(FLG)  GO  TO  50 
A(K)=A(K)+P(J) 

GO  TO  200 

50  L  =  JDIAGdO  "  K 
DO  100  1=1, NEL 
M  =  LD(I) 

IF(M.GT.K  .OR.  MoEO.O)  GO  TO  100 
M  =  L  M 
A(M)=A(M)+S(I, J) 

100  CONTINUE 
200  CONTINUE 
RETURN 
END 
C 

FUNCTION  D0T(A,B,N) 

UECTOR  DOT  PRODUCT 
DIMENSION  A(1),B(1) 

DOT  =  0.0 
DO  100  I=1pN 

100  DOT  =  DOT  +  A(I)^^*B(I) 

RETURN 

END 

C 

SUBROUTINE  PLOAD( ID, F, B, NN, P) 

load  UECTOR  IN  COMPACT  FORM 
DIMENSION  ID(l),F(l),B(l) 

DO  IGO  N=1,NN 
J=ID(N) 

100  IF(J.GT.O)  B(J)=F(N)---P 
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RETURN 

END 

C 

FUNCTION  PROPLD(T»J) 

PROPORTIONAL  LOAD  TABLE  (ONE  LOAD  CARD  ONLY) 

COMMON  /CTDATA/  OjHEAD(20)» NUMNPsNUMELsLAYERjNEQj IPR 
DIMENSION  A(5) 

IF  (J  .LE.  0)  GO  TO  200 
C....  INPUT  TABLE  OF  PROPORTIONAL  LOADS 
1=1 

READ(5»1000)  KsLjTtilNjTMAXc (A(KKK).KKK=1»5) 

L!RITE(G»2000)  OjHEAD^  I»l<f.LpTMIN!iTMA><9  (A(I<KK)»KKK-1>5) 

RETURN 

C....  COMPUTE  UALUE  AT  TIME  T 
200  PROPLD  =0.0 

IF(T.LT.TMIN  .OR.  T.GT.TMAX)  RETURN 
L  =  MAXOCL.l)  , 

PROPLD  =  A(l)-i-A(2)-"T+A(3)^f(SIN(A(4)«T+A(5)))“~L 
RETURN 

1000  FORMAT(2I5,7F10.0)  ^  ,,,,  Kn.Mr,r-D 

2000  FORMAT(Al920A4//5X923HPROPORTIONAL  LOAD  TABLE//1 IN  NUMBER  , 

1  43H  TYPE  EXP.  MINIMUM  TIME  MAXIMUM  TIMEs ISX? BrIAl. 13X> 

2  2HA2p 13X9  2HA39  13X9  2HA49 ISX. 2HA5/(3I89  7G15. 5) ) 

END 

C 

SUBROUTINE  PRTDISCUL. ID9 X9 B9 Ft  T9 NDM9 NDF) 

C«-::-:;-“-  OUTPUT  NODAL  UALUES 
LOGICAL  PCOMP 
COMMON  /PROLOD/  PROP 

COMMON  /CTDATA/  O9 HEAD(20) 9 NUMNPt NUMELt LAYERt NEOt IPR 
COMMON  /LABELS/  PDIS(G) 9 A(G) 9 BC(2) 9 DI (G) 9 CD(3) . FD(3) 

COMMON  /TMDATA/  TIME9 DT9 DDTt FORCE9 ALPHA 

DIMENSION  X(NDM9 1)9B(1)9UL(G)9 IDCNDFt 1)9F(NDF9 l).T(l) 

DATA  BL/4HBLAN/ 

DO  102  N=l9NUMNP 
IF(PC0MP(X(l9N)9BL))  GO  TO  101 
DO  100  I=l9NDF 
UL(I)  =  F(l9N)«PR0P 
K  =  IABS(ID(l9N)) 

100  IF(K.GT.O)  UL(I)=B(K) 

T(N)=UL(3) 

101  CONTINUE 

102  CONTINUE 

WRITEO.EOOl)  (T<I)9l=l9NUMNP) 

RETURN 

2001  F0RMAT(GE12.4) 

END 

SUBROUTINE  FSTREA(UL9XL9LDt  P. IX9 IDtXtF* JDIAG9  DRtBjNDFtNDMtNENt 
NSTtNNEQ) 

ELEfiEHT  ROUT  I  HE 

COMHOH  /CTDATA/  0,  HEAD ( 20 )j.HUHHP9HUMEL»  LAYER* MEQ?  IPR 

COMMOH  /ELDATA/  MpMELs-MCT 

COmOH  /ISWIDX/  ISW 

COmOM  /PROLOD/  PROP  ^ 

DIflEHSIGN  ULCHDF* 1)*XL(HDM9 1)*LD(HDF* 1)*P(1)» IXCNEN* 1)* 

1  IDCNDF*  i)*X(NDM9  l)*F(hDF*  1)*  JDIAGCD* DR(1)*B(1)*S(1) 

IFdSW.EQ.S)  CALL  PLOADC ID* F* DR* HMEQ* PROP) 

MCT=0 

DO  110  H=1*NUHEL 

CALL  PFORn(UL*XL*LD* IX* ID* X*F* B* NDF* NDM*HEH* ISW) 

CALL  ELHT01(UL*XL*IX(1*H)*P*HDF*NDM*NST*ISN) 

IF(ISW„MEc4)  CALL  ADDSTFCDR* S* P*  JDIAG* LD*  1* HEL-^NDF* -FALSEd 
110  COHTIHUE 
RETURH 
END 


SUEROUTIHE  PFORMCUL* XL* LD* IX* ID* X* F* U* HDF* MDH* HEM* ISW) 
FORM  LOCAL  ARRAYS 
COMMOH  /ELDATA/  H*HEL*MCT 
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COMMON  /PROLOD/  PROP 

DIMENSION  UL(NDF»l),XL(NDM»l)»LD(NDFf 1),IX(NEN»1)»ID(NDF»1)» 
^  X(NDM,1),F(NDF,1)pU(1) 

DO  108  I=1»NEN 
II  =  IX(I>N) 

IFCII  .NEo  0)  GO  TO  105 
DO  103  J=1,NDM 

103  XL(J.I)  =  0, 

DO  104  J=l9NDF 
ULCJpI)  =  0. 

104  LD(J.I)  =  0 
CO  TO  103 

105  IID  =  II-NDF  -  NDF 
NEL  =  I 

DO  lOB  J=1.NDM 
lOG  XLCJpI)  =  XCJpII) 

DO  107  J=1,NDF 
K  =  IABS(ID(J9lI)) 

ULCJ^I)  =  FCJ^ID^fPROP 
IF(K.GToO)  UL(J9l)=U(K) 

IF(ISWoEQ.G)  K=IID+J 

107  LDCJ^I)  =  K 

108  CONTINUE 
RETURN 
END 

C 

SUBROUTINE  ELMTOl (UL* XL, IX, P» NDF, NDM^ NST, ISW) 

LINEAR  ELASTIC  IN-PLANE  BENDING  ELEMENT  ROUTINE 
LOGICAL  TAN 

COMMON  /ELDATA/  N,NEL,MCT 

COMMON  /MTDATA/  RHO, UU12, El, E2, G12, G13, GE3, THK, WIDTH 
COMMON  /COMPST/  ABD(B,B),DS(2,2),QBR(3,3,25),QBS(2,2,25)f 
-  TH(25),ZK(25) 

COMMON  /DMATIX/  D( 10) , DB(G, B) , LINT 
COMMON  /TMDATA/  TIME, DT, DDT, FORCE,  ALPHA 
COMMON  /GAUSSP/  SG( IG) , TG( IB) , WG( IG) 

COMMON  /EXTRAS/  TAN 

DIMENSION  UL(NDF,1),XL(NDM, 1),IX(1),P(1),SHP(3,12), 

1  SIGT(3),SIGB(3),SIGS(2),EPT(3),EPB(3),EPS(2) 

C 

DO  20  L=1,NST 
20  P(L)  =  0.0 

C..OO  COMPUTE  NEUTRAL  STRAINS  AND  STRESS  RESULTANTS 
L  =  D(l) 

IF(ISW.EQ.4)  L=D(3) 

CALL  PGAUSS(L,LINT) 

DO  GOO  L=1,LINT 

C  ..  COMPUTE  ELEMENT  SHAPE  FUNCTIONS 

CALL  SHAPE(SG(L) , TG(L) , XL, SHP, XSJ, NDM, NEL, IX, .FALSE, ) 

C  .0  COMPUTE  STRAINS  AND  COORDINATES 
DO  410  1=1,3 
EPT(I)  =  OoO 
410  EPB(I)  =  OoO 
DO  420  1=1,2 
420  EPSd)  =  0.0 
XX  =  OoO 
YY  =  0.0 
DO  430  J=1,NEL 
XX  =  XX  +  SHP(3,  J)-?:-XL(l,  J) 

YY  =  YY  +  SHP(3,  J)^>XL(2,  J) 

C  ..  IN-PLANE  STRAINS 

EPT(l)  =  EPT(l)  +  SHPCl,  J)^:'UL(1,J) 

EPT(2)  =  EPT(2)  +  5HP(2, J)^UL(2, J) 

EPT(3)  =  EPT(3)  +  SHPCl,  J)-UL(2,  J)  +  SHP(2,  J)«-UL(  1,  J) 

C  ..  BENDING  CURUATURES 

EPB(l)  =  EPB(l)  -  SHPCl, J)*»UL(4, J) 

EPB(2)  =  EPB(2)  -  SI:P(2, J)-ULC5, J) 

EPB(3)  =  EPB(3)  “  5HP(1, J)^UL(5, J)  -  SHPC2, J)*UL(4, J) 

C  ..  SHEARING  STRAINS 

EPSd)  =  EPSd)  +  SHPCl, J)*UL(3, J)  -  SHP(3, J)*UL(4, J) 


PFOR  4 
PFOR  5 
PFOR  B 
PFOR  7 
PFOR  8 
PFOR  3 
PFOR  10 
PFOR  11 
PFOR  12 
PFOR  13 
PFOR  14 
PFOR  15 
PFOR  IG 
PFOR  17 
PFOR  18 
PFOR  19 
PFOR  20 
PFOR  21 
PFOR  22 
PFOR  23 
PFOR  24 
PFOR  25 
PFOR  2G 
PFOR  27 
PFOR  28 
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ELMT  23 
ELMT  24 
ELMT  25 
ELMT  2S 
ELMT  27 
ELMT  28 
ELMT  29 
ELMT  30 
ELMT  31 
ELMT  32 
ELMT  33 
ELMT  34 
ELMT  35 
ELMT  3G 
ELMT  37 
ELMT  38 
ELMT  39 
ELMT  40 
ELMT  41 
ELMT  42 
ELMT  43 
ELMT  44 
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430  EPSO)  =  EPS(2)  +  SHP(2» J)»UL(3, J)  -  SHP(3. J)*UL(5, J) 
IFCISW.EQ.S.AND.TAN) 

WRITE(9,9001)  M»Lp(EPB(II), 11=1.3), (EPS(II). 1=1.2) 

9001  F0RmT(2IG,5E12.4) 

C  ..  COMPUTE  3TRESS  RESULTANTS 
DO  440  1=1.3 
SIGT(I)  =  0. 

SIGB(I)  =  0. 

DO  440  J=1.3 

SIGT(I)  =  SIGT(I)  +  ABDC!. J)«EPT(J)  +  ABD( I. J+3)*EPB(J) 

440  SIGB(I)  =  SIGB(I)  +  ABD(  1+3.  J)-::EPT( J)  +  ABD(  1+3, J+3)«-EPB( J) 
DO  450  1=1,2 
SIGS(I)  =  0. 

DO  450  J=1.2 

450  SIGS(I)  =  SIGSCI)  +  DSd. J)»EPS(J) 

IF(ISW.GT.4)  GO  TO  B20 

C  ..  OUTPUT  STRESS  RESULTANTS  AND  STRAINS 
MCT  =  MCT  -  2 
IF(MCT.GT.O)  GO  TO  470 
WRITE(G.2001)  TINE 


470 


MCT  =  50 

WRITE(B,2002)  N.XX.YY.EPT. 


EPB. EPS. SIGT, SIGB, SIGS 


GO  TO  GOO 

C....  COMPUTE  INTERAL  FORCES 
GEO  DU  =  XSJ«LIG(L) 

J1  =  1 


DO  GIO  J=1.NEL 
P(J1  )  =  F(J1  ) 
P(J1+1)  =  P(J1+1) 
P(Jl+2)  =  P(Jl+2) 
P(Jl+3)  =  P(Jl+3) 

P(Jl+4)  =  P(Jl+4) 


-  (SHPCl, J)«SIGT(1)+SHP(2,J)«SIGT(3))»DU 

-  (SHP(2. J)^=^SIGT(2)+SHP(1, J)«-SIGT(3))»DU 

-  (SHPd. J)»SIGS(1)+SHP(2, J)«SIGS(2))«DU 

+  (SHPCl.  J)«SIGB(1)+SHP(2.  J)«-SIGB(3)+SHP(3,J) 
^JSIGSd))-DU 

+  (SHP(2, J)»SIGB(2)+SHPd. J)«SIGB(3)+SHP(3,J) 
»SIGS(E))*DU 


GIO  J1  =  J1  +  NDF 
GOO  CONTINUE 
RETURN 

C 


2001  FORMAT! IHl// 

^  5X.GHTIME  =.E12.3//5X,33HELEMENT  STRAINS/STRESS  RESULTANTS// 

1  8H  ELEMENT. 3X.7H1-C00RD.3X.7H2-C00RD,4X,9HXX-STRAIN,4X, 

2  9HYY-STRAIN. 4X. 9HXY-STRAIN, 3X, lOHKXX-STRAIN, 3X. 

3  10HKYY-STRAIN,3X,10HKXY-STRAIN,4X,9HSX-STRAIN.4X. 

4  9HSY-STRAIN/28X. 8(GX. 7H-STRESS)/) 

2002  FDRMATd8.2F10.4,8E13.4/28X.8E13.4) 

END 


C 

SUBROUTINE  PGAUSS(LL,LINT) 

GAUSSIAN  POINTS  AND  WEIGHTS  FOR  TWO  DIMENSIONS 
COMMON  /GAUSSP/  SGCIG) . TG( IG) , WGCIG) 

DIMENSION  LR(S).LZ(3),LW(9).WR(2),GR(E).GC(2) 

DATA  LR/-1. l.l.-l. 0. 1,0.-1.0/,L2/-1,-1. 1. 1. -1. 0. 1. 0. 0/ 
DATA  LW/4-“25,4'-f40,S4/ 

DATA  GR/0 . 8G113G31 1594053, 0 . 33398104358485G/ 

DATA  GC/ 1.0. 0.3333333333/ 

DATA  WR/0 . 347854845137454, 0 . G521451548G254B/ 

LINT  =  LL«LL 
L=IABS(LL) 

GO  TO  (1.2.3,4).L 
C....  1X1  INTEGRATION 

1  SG(l)  =  0. 

TGCl)  =  0. 

WGd)  =  4. 

RETURN 

C....  2X2  INTEGRATION 

2  G  =  l./SQRT(3.) 

IF(LL.LT.O)  G=l. 

DO  21  1=1.4 
SGd)  =  G*LRd) 

TGCI)  =  G«LZd) 


ELMT 

45 

ELHT 

4B 

ELni 

47 

ELliT 

4B 

ELNT 

49 

ELNT 

SO 

ELNT 

Si 

ELNT 

£•; 

L-. 

ELNT 

ELNT 

54 

ELNT 

55 

ELNT 

r':,p 

ELNT 

s  / 

ELNT 

SB 

ELNT 

59 

ELNT 

GO 

ELNT 

G1 

ELNT 

G9 

ELNT 

G3 

ELNT 

G4 

ELNT 

G5 

ELNT 

GG 

ELNT 

G7 

ELNT 

Go 

ELNT 

G9 

ELNT 

70 

ELNT 

71 

ELNT 

7L 

ELNT 

73 

ELNT 

74 

ELNT 

75 

ELNT 

7B 

ELNT 

77 

ELNT 

70 

ELNT 

79 

ELNT 

SO 

ELNT 

81 

ELNT 

B2 

ELNT 

83 

ELNT 

84 

ELNT 

85 

ELNT 

8S 

ELNT 

87 

ELNT 

88 

ELNT 

83 

ELNT 

90 

PGAU 

i 

PGAU 

2 

PGAU 

o 

PGAU 

4 

PGAU 

5 

PGAU 

G 

PGAU 

7 

PGAU 

8 

PGAU 

3 

PGAU 

10 

PGAU 

11 

PGAU 

12 

PGAU 

13 

PGAU 

14 

PGAU 

15 

PGAU 

IG 

PGAU 

17 

PGAU 

18 

PGAU 

13 

PGAU 

20 

PGAU 

21 

PGAU 

22 

PGAU 

23 

137 


21  WG(I)  =  1. 

RETURN 

C,...  3X3  INTEGRATION 

3  G  =  SQRT(0,6) 

IF(LL.LT.O)  G=l. 

H  =  lo/81. 

DO  31  1=1,9 
SG(I)  =  G-J'-LRCI) 

TG(I)  =  G^'LZd) 

31  UG(I)  =  H-^LW(I) 

RETURN 

C....  4X4  INTEGRATION 

4  DO  41  1=1,4 

11  =  1+M0D(I+1,2) 

12  =  1 

IF(I.GTo2)  12  =  2 
DO  41  J=l,4 
JJ  =  (I-l)«4+J 
SG(JJ)  =  LR(J)^fGR(Il) 

IF(LL.LT.O)  SG(JJ)  =  LRCJ)*GC(I1) 

TG(JJ)  =  LZ(J)«*GR(I2) 

IF(LLoLT.O)  TG(JJ)  =  LZ(J)*GC(I2) 

41  WG(JJ)  =  NR(I1)«WR(I2) 

RETURN 

END 

C 

SUBROUTINE  SHAPECSS, TT, X,SHP,X5J,NDM,NEL, IX,FLG) 

SHAPE  FUNCTION  ROUTINE  FOR  TWO  DIMENSIONAL  ELEMENTS 
LOGICAL  FLG 

DIMENSION  SHP(3,4),X(NDM, 1),S(4),T(4),XS(2,2),SX(2,2),IX(9) 
DATA  S/-0,5,0o5,0.5,-0.5/,T/-0.5,-0.5,0.5,0.5/ 

C.P..  FORM  4--N0DE  QUADRILATERIAL  SHAPE  FUNCTIONS 
DO  100  1=1,4 

SHPO,!)  =  (0o5+S(I)«SS)*(0.5+T(I)^TT) 

SHP(1,I)  =  S(I)^K0o5+T(I)^^TT) 

100  SHP(2,I)  =  T(I)*(0.5+S(I)«-S5) 

IF(NELoGE,4)  GO  TO  120 

C _  FORM  TRIANGLE  BY  ADDING  THIRD  AND  FOURTH  TOGETHER 

DO  110  1=1,3 

110  SHP(I,3)  =  SHP(I,3)+SHP(I,4) 

C.*..  ADD  QUADRATIC  TERMS  IF  NECESSARY 

120  IF(NELcGTc4  ,ANDo  NEL.LT.IO)  CALL  SHAP2(SS, TT, SHP, IX, NED 
C..O-  ADD  CUBIC  TERMS  IF  NECESSARY 

IF(NEL.GTc9)  CALL  SHAP3(SS, TT, SHP, IX, NED 
C...,  CONSTRUCT  JACOBIAN  AND  ITS  INUERSE 
DO  130  1=1, NDM 
DO  130  J=l,2 
XS(I,J)  =  OoO 
DO  130  K=1,NEL 

130  XS(I,J)  =  XS(I,J)+  X(J,K)*SHP(I,K) 

XSJ  =  X3(l,  l)wXS(2,2)“XS(l,2)»‘XS(2,l) 

IFCXSJ  .GT.  OoOOOOOOOl)  GO  TO  135 
WRXTECGpEOOO)  IX 
C  STOP 

135  IF(FLG)  RETURN 

SX(i,l)  =  XS(E,2)/XSJ 
SX(2,2)  =  X5(1,1)/XSJ 
SX(1p2)  =  -XS(1,2)/XSJ 
SX(2,1)  =  ~XSCE,1)/XSJ 
Coo*.  FORM  GLOBAL  DERIUATIUES 
DO  140  1=1, NEL 

TP  =  5HP(1,I)‘-SX(1,1)+SHP(2,I)-'SX(2,1) 

SHP(E,I)  =  SHP(1,I)^-SX(1,2)+SHP(2,I)-^SX(2,2) 

140  S!-IP(1,I)  =  TP 
RETURN 

2000  FORMAT (5X,67H^‘-”SHAPE  ERROR  01--  ZERO  OR  NEGATIUE  JACOBIAN  DET. 
-ELEMENT  NODESs /20X, 1214) 

END 

C 

SUBROUTINE  SHAP2(S, T, SHP, IX, NED 
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100 

C  t  •  •  • 


101 

102 

103 

C  •  •  •  • 

104 

C  •  *  •  • 
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100 
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103 


ADD  QUADRATIC  FUNCTIONS  AS  NECESSARY 
DIMENSION  IX(9).SHP(3»12) 

S2  =  (l.-S»S)/2. 

T2  =  (l.-T-»T)/20 
DO  100  1=5. NEL 
DO  100  J=1.3 
SHPCJsI)  =  0.0 

MIDSIDE  NODES  (SERENDIPITY) 
IF(IX(5).EQ.O)  GO  TO  101 
SHPd.S)  =  -S«(l.-T) 

SHP(2,5)  =  -S2 

SHPO.S)  =  S2*(l.-T) 

IF(NEL.LT.G)  GO  TO  107 
IF(IX(G).EQ.0)  GO  TO  102 
SHPCl.G)  =  T2 
SHPCEjG)  =  -PK1.+S) 

SHPO.G)  =  T2<J(1.+S) 

IF(NEL.LT,7)  GO  TO  107 
IF(IX(7).EQ.0)  GO  TO  103 
SHP(l07)  =  -S«(l.+T) 


SHP(297)  =  S2 
SI-IP(3,7)  =  S2-;K1.+T) 

IFCNEL.LT.S)  GO  TO  107 
IF(IX(8).EQ.O)  GO  TO  104 
SHPd.S)  =  -T2 
SHP(2.8)  =  -TK-d.-S) 

SHP(3.8)  =  T2»d.-S) 

INTERIOR  NODE  (LAGRANGIAN) 

IF(NEL.LT.3)  GO  TO  107 
IF(IX(S).EQ.O)  GO  TO  107 
SHPd.9)  =  -4.«S«T2 
SHP(2.9)  =  -4.»T«S2 
SHP(3.9)  =  4.*S2*T2 

CORRECT  EDGE  NODES  FOR  INTERIOR  NODE(LAGRANGIAN) 

DO  lOG  J=1.3 
DO  105  1=1.4 

SHP(J.I)  =  SHP(J.  I)  -  0.25:fSHP(J.9) 

DO  lOG  1=5.8 

IFdX(I).NE.O)  SHP(J.I)  =  SHP(J.I)  -0.5»SHP( J. 9) 

CORRECT  CORNER  NODES  FOR  PRESENCE  OF  MIDSIDE  NODES 


K  “  8 

DO  103  1=1.4 
L  =  I  +  4 
DO  108  J=1.3 

SHP(J.I)  =  SHPCJ.I)  -  0.5»(SHP(J.K)+SHP(J.L)) 

K  =  L 

RETURN 

END 


SUBROUTINE  SHAP3(S. T. SHP. IX. NED 
ADD  CUBIC  FUNCTION  AS  NECESSARY 
DIMENSION  IX(1E)pSHP(3. 12) 

DO  100  1=5. NEL 
DO  100  J=1,3 
SHP(J.I)=0.0 
IF(IX(5)-EQ*0)  GO  TO  101 
Sl=“l./3. 

Tl=“l. 

CALL  CSHAPECS.T.Sl.Tl.SHP.l.S) 
IF(IX(6).EQ,0)  GO  TO  102 
Sl=l. 

Tl=-l./3o 

CALL  CSHAPE(S.T.S1.T1pSHP.2.6) 
IF(IX(7)oEQ.0)  GO  TO  103 
Sl=l./3. 

Tl  =  i. 

CALL  CSHAPECSpT.SIpTIpSHP.I.Z) 
IF(IX(8)oEQoO)  GO  TO  104 
Sl=-1, 

Tl=lo/3. 


(SERENDIPITY) 
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CALL  CSHAPE(S,T,S1,T1,SHP»2,8) 

104  IF(IX(9).EQ.O)  GO  TO  105 

Tl=-lo/3« 

CALL  CSHAPE(S»T.S1,T1pSHP^2j9) 

105  IFCNELoLTolO)  GO  TO  EOO 
IF(IX(10).EQ.0)  GO  TO  lOG 
91=1. /3. 

Tl=“l. 

CALL  CSHAPE(S9T,SlfTl»SHP, 1,10) 
lOG  IF(hEL.LT.ll)  GO  TO  200 
IFdXClD.EQoO)  GO  TO  107 
91=1. 

Tl=l./3. 

CALL  C9HAPE(5,T,Sl,Tl93HP,2,ll) 

107  IF(NEL.LT.12)  GO  TO  200 
IF(IX(12).EQ.O).  GO  TO  200 
91=~l./3. 

Tl  =  l. 

CALL  C5HAPE(99T,5l9Tl,9HP,l,12) 

C.i..  CORRECT  CORNER  NODES 
200  DO  210  1=1,4 
11=1+4 
12=1+8 

IF(I.EQ.l)  13=1+7 
IF(I.GT.l)  13=1+3 
IF(I.LT.4)  14=1+9 
IF(I.EQ.4)  14=1+5 
DO  210  J=l,3 

210  3HP(J,I)=5HP(J, I)“2./3.*(SHP(J, I1)+SHP(J,I2))-1./3.*(SHP(J,  13) 
-  +SHP(J,I4)) 

RETURN 

END 

C 

SUBROUTINE  CSHAPECS, T, SI, T1 , SHP, K, L) 

SUPPLEMENTAL  ROUTINE  FOR  THE  SHAPE  FUNCTIONS 
DIMENSION  SHP(3,12) 

C=9o/32. 

GO  TO  (1,2),K 

1  SHP  ( 1 ,  L )  =C*  ( 1 .  +T 1  *T )  *  ( 9 .  *S  1  -2 .  *S-27 .  **^91*9*9 ) 

SHP  ( 2 ,  L )  =C-''^T  1  ( 1 . -S*-'S )  ( 1 . +9 .  S 1  ^^S ) 
SHP(39L)=C^(lo+Tl*T)«-(l.*-S«S)*(l.+9.*Sl*S) 

RETURN 

2  SHP(l,L)=C-»SP^(l.-Ti^'T)^^(l.+9.*Tl^T) 
SHP(2,L)=C-5^(l.+Sl^^S)-(9.^Tl-2.^tT’-27.^Tl«T*T) 
5HP(3,L)=C*(l.+Sl*S)’»(l.“T*T)^(l.+9.-»Tl*T) 

RETURN 

END 

C 

SUBROUTINE  PMESH(IDL,XL, IX, ID, X, F, JDIAG,NDF,NDM,NEN,NKM) 
INPUT  MESH  DATA 
LOGICAL  PRT,ERR,PCOMP 

COMMON  /CTDATA/  0, HEAD (20) , NUMNP, NUMEL, LAYER, NEQ, IPR 
COMMON  /MTDATA/  RH0,UU12, El, E2,G12,G13,G23,THK, WIDTH 
COMMON  /LABELS/  PDIS(G) , A(B) , BC(2) , DI (G) , CD(3) , FD(3) 

COMMON  /EXDATA/  QLAW(4) 

COMMON  /RODATA/  UR,IQ,ND3 

DIMENSION  IDL(G),XL(7), IXCNEN, 1), IDCNDF, 1),X(NDN, 1), 

^  F(NDF9l),DUM(l),WD(13), JDIAG(l) 

DATA  WD/4HC00R,4HELEM,4HMATE,4HB0UN,4HF0RC,4HR0D  , 

^  4HEND  9  4HPRIN, 4HN0PR, 4HPAGE, 4HEXPE/ 

DATA  BL/4HBLAN/, LIST/11/, PRT/ . TRUE . / 

C....  INITIALIZE  ARRAYS 

ERR  =  .FALSE. 

DO  501  1=1,4 
501  QLAW(I)=0. 

DO  502  N=1,NUMNP 
DO  502  1=1, NDF 
ID(I,N)=0 
F(I,N)=0. 
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SHAP  50 
SHAP  51 
SHAP  52 
SHAP  53 
SHAP  54 

CSHA  1 
CSHA  2 
CSHA  3 
CSHA  4 
CSHA  5 
CSHA  G 
CSHA  7 
CSHA  8 
CSHA  9 
CSHA  10 
CSHA  11 
CSHA  12 
CSHA  13 
CSHA  14 

PMES  1 
PMES  2 
PMES  3 
PMES  4 
PMES  5 
PMES  G 
PMES  7 
PMES  8 
PMES  9 
PMES  10 
PMES  11 
PMES  12 
PMES  13 
PMES  14 
PMES  15 
PMES  16 
PMES  17 
PMES  18 
PMES  19 
PMES  20 
PMES  21 


n  n 


502  CONTINUE 

C _  READ  A  CARD  AND  COMPARE  WITH  MACRO  LIST 

10  READ(5,1000)  CC 
DO  20  I=l5LIST 

20  IF(PC0MP(CC4JD(I)))  GO  TO  30 
GO  TO  10 

30  GO  TO  (1<.2p3p4»5.Gj?.8jSi.  II9 12).  I 
C....  NODAL.  COORDINATE  DATA  INPUT 

1  DO  102  N=1.NUMNP 
102  X(1,N)=  BL 

CALL  GENUEC(NDM.XL.X.CD.PRT.ERR) 

GO  TO  10 

C....  ELEMENT  DATA  INPUT 

2  L=0 

DO  20G  I=1.NUMEL.50 

IF(PRT)  WRITE(E.2001)  O.KEAD.  (K.K=1.NEN) 

J  =  MIN0(NUMEL.I+49) 

DO  EOb  N=I.J 
IF(L-N)  200.202.203 

200  READ(5.1001)  L. (IDL(K).K=1.NEN).LX 
IF(L.EQoO)  L=NUMEL+i 
IF(LX.EQ.O)  LX=1 

IF(L-N)  201.202.203 

201  WRITE(6.3001)  L.N 
ERR  =  .TRUE. 

GO  TO  EOG 

202  NX  =  LX 

DO  207  K=1.NEN 
207  IX(K,L)  =  IDL(K) 

GO  TO  205 

203  IXCNEN.N)  =  IX(NEN.N-l) 

DO  204  K=1.NEN 

IX(K.N)  =  IX(K.N-l)  +  NX 

204  IF(IX(K.N-l).EQ.O)  IX(K.N)  =  0 

205  IF(PRT)  WRITE(G.EOOE)  N. (IX(K.N).K=1.NEN) 

EOG  CONTINUE 

GO  TO  10 

C....  MATERIAL  DATA  INPUT 

3  WRITE(G.2004)  O.HEAD 
CALL  NATL IB 

GO  TO  10 

C....  READ  IN  THE  RESTRAINT  CONDITIONS  FOR  EACH  NODE 

4  IF(PRT)  WRITECG.EOOO)  O.HEAD. (I. BC. I=1»NDF) 

N  =  0 

NG  =  0 
420  L  =  N 
LG  =  NG 

READ(5.1001)  N.NG.IDL 

IFCN.LE.O  .OR.  N.GT.NUMNP)  GO  TO  50 

DO  41  1=1. NDF 

41  IFCl'.nLo  .and.  IDLCD.EQ.O  .and.  ID(I.L).LT.O)  ID(ItN)=-l 
LG  =  ISIGN(LG.N-L) 

42  L  =  L+LG 

IF((N-L)’>LG  .LE.  0)  GO  TO  420 
DO  43  1=1, NDF 

43  IF(ID(I.L-LG)  .LT.  0)  ID(I.L)  =  -1 
GO  TO  42 

50  DO  48  N=1.NUMNP 
DO  4G  1=1. NDF 

4G  IFdDCI.N)  .NE.  0)  GO  TO  47 
GO  TO  48 

47  IF(PRT)  WRITE(B.2007)  N. (ID( I. N) . 1=1, NDF) 

48  CONTINUE 
GO  TO  10 

C....  FORCE/DISPL  DATA  INPUT 

5  CALL  GENUECCNDF.XL.F.FD.PRT.ERR) 

GO  TO  10 

END  OF  MESH  DATA  INPUT 
COMPUTE  THE  PROFILE  OF  GLOBLE  ARRAYS 


PMES  22 
PMES  23 
PI1ES  24 
PMES  2b 
PMES  23 
PMES  27 
PMES  £8 
PiiE:3  2S 
PMES  30 
Pi1ES  31 
PMES  32 
PMES  33 
PMES  34 
PMES  35 
FiTES  3S 
FMES  37 
PMES  33 
FI1ES  Li'J 
PMES  40 
PMES  41 
PMES  4E 
PMES  4-3 
PMES  44 
PiIlS  45 
PMES  4B 
PMES  47 
PMES  48 
PMES  49 
PMES  50 
PMES  5i 
PMES 
PMES  53 
PMES  54 
PMES  55 
PMES  5S 
PMES  57 
PMES  53 
PMES  59 
PMES  GO 
PMES  G1 
PMES  G2 
PMES  63 
PMES  G4 
PMES  G5 
PMES  GG 
PMES  67 
PMES  63 
PMES  B3 
PMES  70 
PMES  71 
PMES  72 
PMES  73 
PMES  74 
PMES  75 
PMES  7G 
PMES  77 
PMES  78 
PMES  79 
PMES  80 
PMES  81 
PMES  82 
PMES  03 
PMES  84 
PMES  85 
PMES  09 
PMES  87 
PMES  88 
PMES  89 
PMES  SO 
PMES  91 
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7  IF (ERR)  STOP 

CALL  PROFIL(JDIAG,ID> IX»NDF»NENjNKMf PRT) 

RETURN 

C.o.o  PRINT  OPTION 

8  PRT  =  .TRUE. 

GO  TO  10 

C....  NOPRINT  OPTION 
3  PRT  =  .FALSE. 

GO  TO  10 

C....  READ  IN  PAPER  EJECTION  OPTION 

11  READ(5,1000)  0 
GO  TO  10 

C....  INPUT  EXPERIMENTAL  INDENTATION  LAW 

12  READ(5,1007)  (QLAW( I ) t I=l» 4) 

WRITE(B,2008)  OtHEAD? (QLAW( I) > I=l» 4) 

GO  TO  10 

C....  INPUT  INITIAL  IMPACT  CONDITION 
G  WRITE(6,2003)  O^HEAD 
READ(5,1002)  NO, INDFpUR 
WRITE(G.2010)  NQtINDFfUR 
F(INDF,NQ)=1.0 
IQ=ID(INDF.NQ) 

GO  TO  10 

C _  INPUT/OUTPUT  FORMATS 

1000  F0RMAT(A4,75X,A1) 

1001  FORMATCIGIS) 

1002  FORMAT(2I59F10.0) 

2000  FORMAT(A1i>26a4//5X»  lOHNODAL  B.C. » 7X//BXf  SHNODE^  ,9(I7f  A4» A2)/1X) 

2001  FORMATC A1 , 20A4//5X9  8HELEMENTS//3Xj  7HELEMENT » 

^  14(13. 5H  NODE)/(20X. 14(13. 5H  NODE))) 

2002  FORMAT(I10.14I8/(10X.14I8)) 

2004  FORMAT(A1.20A4//5X. 13HMATERIAL  PROPERTIES) 

2007  FORMAT(I10.3I13) 

2008  FORMAT(A1,20A4//5X,?^EXPERIMENTAL  INDENTATION  LAW?^// 

1  10X.?^CONTACT  COEFFICIENTS  ?i.E12.4/ 

2  10X.?^CRITICAL  INDENTATION:  E12.4/ 

3  10X.?^CONSTANT  Ss  E12.4/ 

4  10X.?^POWER  INDEX  OF  UNLOADING  LAW:^  F12.3) 

3001  FORMAT (5X.26H-‘«PMESH  ERROR  01^^*^^  ELEMENT.  15. 

-  22H  APPEARS  AFTER  ELEMENT. 15) 

2003  F0RMAT(A1.20A4.//5X.?^IMPACT  OF  LAMINATED  PLATEj^) 

2010  FORMAT(//10X.?^IMPACT  NODAL  POINT:  j^.IlO/ 

-  10X.?^IMPACT  D.OoFoS  ?^.I10/ 

-  lOX.^INITIAL  IMPACT  UELOCITY: E12.4) 

END 

C 

SUBROUTINE  GENUEC(NDM. XL. X. CD. PRT. ERR) 

GENERATE  REAL  DATA  ARRAYS  BY  LINEAR  INTERPOLATION 
LOGICAL  PRT.ERR.PCOMP 

COMMON  /CTDATA/  0. HEAD(20) . NUMNP. NUMEL. LAYER. NEQ. IPR 
DIMENSION  X(NDM. 1).XL(7).CD(3) 

DATA  BL/4HBLAN/ 

N=0 

NG=0 

102  L=N 
LG=NG 

READ(5.1000)  N.NG.XL 

IFCN.LE.O  .OR.  N.GT.NUMNP)  GO  TO  108 

DO  103  1=1. NDM 

103  X(I.N)=XL(I) 

IF(LG)  104.102.104 

104  LG=ISIGN(LG.N-L) 

LI =( I ABS ( N-L+LG ) -1 ) /I ABS (LG ) 

DO  105  1=1. NDM 

105  XL(I)=(X(I.N)-X(I.L))/LI 
iOG  L=L+LG 

IF((N~L)^'LG  .LE.  0)  GO  TO  102 
IF(L.LE.0  .OR.  L.GT.NUMNP)  GO  TO  110 
DO  107  1=1. NDM 


PMES  92 
PMES  93 
PMES  94 
PMES  95 
PMES  3S 
PMES  97 
PMES  98 
PMES  99 
PMESlOO 
PMESlOl 
PMES102 
PMES103 
PMES104 
PMES105 
PMESlOB 
PMES107 
PMESIOS 
PMES109 
PMESllO 
PMESlll 
PMES112 
PMES113 
PMES114 
PMES115 
PMES116 
PMESil7 
PMES118 
PMES119 
PMES120 
PMES121 
PMES 122 
PMES 123 
PMES124 
PMES 125 
PMES12G 
PMES127 
PMES128 
PMES 129 
PMES130 
PMES131 
PMES 132 
PMES133 
PMES134 
PMES135 
PMES13G 
PMES 137 

GENU  1 
GENU  2 
GENU  3 
GENU  4 
GENU  5 
GENU  G 
GENU  7 
GENU  8 
GENU  9 
GENU  10 
GENU  11 
GENU  12 
GENU  13 
GENU  14 
GENU  15 
GENU  IB 
GENU  17 
GENU  18 
GENU  19 
GENU  20 
GENU  21 
GENU  22 
GENU  23 
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107  X(I*L)=XCI,L-LG)+XL(I) 

GO  TO  lOB 

110  WRITE(B>3000)  L» (CD(I)» I=1p3) 

ERR  =  .TRUE. 

GO  TO  102 

108  DO  103  I=UMUmPp50  .  ^ 

IF(PRT)  WRITE (09 2000 )O^HEAD> (CD(L)9L=lf 3)> (L9CD(l)9CD(2)9L=l9MDri) 

M  =  MINOChUmPpIMB) 

DO  103  J=l9H 

IFCPC0MP(X(l9 J)9BL)  .AMD.  PRT)  URITECB? 2008)  N  . 

109  IFC.HOT.PCOMPCXCIp J)9BL).AriD.PRT)  WRITE(692003)  J9 (XCL? J) fL-l? NDN) 
RETURM 

1000  FORnAT(2I5p7F10.0)  _ _ _ 

2000  FORMAT(Al920A4//5X9  SHMODALp 3A4//6X9 4HN0DE9 9( 17? A49 A2) ) 

2008  FORflATCSXpPlH-^-^CENUEC  WARNING  OI'^^pIIOp 
^  32H  HAS  NOT  BEEN  INPUT  OR  GENERATED) 

3000  FORMAT (5X?44H^--GENUEC  ERROR  Ol^^^ATTEMPT  TO  GENETATE  NODE9  ISp 
1  3H  IN93A4) 

END 


SUBROUTINE  PROFIL( JDIAG, IDp IXp NDFp NENp NKMp PRT) 
COMPUTE  PROFILE  OF  GLOBAL  ARRAYS 
LOGICAL  PRT 

COMMON  /CTDATA/  0, HEAD(20) 9 NUMNP, NUMELp LAYER9NEQ9 IPR 
DIMENSION  JDIAGCDp  IDCNDFp  1)9  IXCNENp  l)f  EQ(2) 

DATA  EQ/4H  DOFpEH.  / 

C....  set  UP  THE  EQUATION  NUMBERS 
NEQ  =  0 

DO  50  N=l9NUMNP 
DO  40  I=1pNDF 
J  =  IDCIpN) 

IF(J)  30920,30 
20  NEQ  =  NEQ  +  I 
ID(I,N)  =  NEQ 
JDIAG(NEQ)  =  0 
GO  TO  40 
30  ID(I,N)  =  0 
40  CONTINUE 
50  CONTINUE 

IF(oNOT.PRT)  GO  TO  70 
WRITE(6,2000)  OpHEAD, (I, EQ, I=1» NDF) 

DO  BO  I=1,NUMNP 

GO  WRITE(G,2001)  I, (ID(K, DpK^IpNDF) 

C....  COMPUTE  COLUMN  HEIGHTS 
70  DO  500  N^IpNUMEL 
DO  400  1=1, NEN 
II  =  IX(I,N) 

IF(II  .EQ.  0)  GO  TO  400 
DO  300  K=1,NDF 
KK  =  IDCKpII) 

IF(KK.EQ.O)  GO  TO  300 
DO  200  J=I,NEN 
JJ  =  IX(J,N) 

IF(JJ.EQ.O)  GO  TO  200 
DO  100  L=1,NDF 
LL  =  ID(L,JJ) 

IF(LL.EQ.O)  GO  TO  100 
M  =  MAX0(KK,LL) 

JDIAG(M)  =  MAXO(JDIAG(M),IABS(KK*"LL)) 

100  CONTINUE 
200  CONTINUE 
300  CONTINUE 
400  CONTINUE 
500  CONTINUE 

C....  COMPUTE  DIAGONAL  POINTERS  FOR  PROFILE 
NKM  =  1 
JDIAG(l)  =  1 
IF(NEQ.EQ.l)  RETURN 
DO  GOO  N=29NEQ 


GENU  24 
GENU  25 
GENU  2B 
GENU  2?' 
GENU  2B 
GENU  20 
GENU  30 
GENU  31 
GENU  32 
GENU  33 
GLr-;u  3^- 
GENU  35 
GENU  36 
GENU  37 
GENU  38 
GENU  39 
GENU  40 
GENU  41 
GENU  42 
GENU  43 

PROF  1 
PROF  2 
PROF  3 
PROF  4 
PROF  5 
PROF  G 
PROF  ? 
PROF  8 
PROF  9 
PROF  10 
PROF  il 
PROF  12 
PROF  13 
PROF  14 
PROF  15 
PROF  IG 
PROF  17 
PROF  18 
PROF  19 
PROF  20 
PROF  21 
PROF  22 
PROF  23 
PROF  24 
PROF  25 
PROF  2B 
PROF  27 
PROF  28 
PROF  23 
PROF  30 
PROF  31 
PROF  32 
PROF  33 
PROF  34 
PROF  35 
PROF  3B 
PROF  37 
PROF  38 
PROF  3^3 
PROF  40 
PROF  41 
PROF  42 
PROF  43 
PROF  44 
PROF  45 
PROF  4B 
PROF  47 
PROF  40 
PROF  49 
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GOO  JDIAG(N)  =  JDIAG(N)  +  JDIAG(N“1)  +  1 
MKM  =  JDIAG(NEQ) 

SOOO  FORMATCAlf 20A4//5X» 16HEQUATI0N  NUMBERS//BX»5HN0DE  » 

-  9(I5,A4,A2)/1X) 

2001  FORMATCIlO^SIll) 

RETURN 

END 

C 

SUBROUTINE  MATLIB 

MATERIAL  PROPERTIES  ROUTINE 
COmON  /CTDATA/  0,  HEAD(20)»NUNNP»NU[1EL9LAYERf  NEQf  IPR 
COMMON  /MTDATA/  RHO, UU12. El fE2»G12,G13,G23,THK, WIDTH 
COMMON  /COMPST/  ABD(G» G) 9 DS(2? 2) » QBR(3p 3f 25) > QBS(29 2»25)» 

-  TH(25)9ZK(25) 

COMMON  /DMATIX/  D(10)9DB(G9G)9LINT 
DIMENSION  WD(5) 

DATA  WD/GH  ISO-9 GH  □RTH09 6HTROPIC9 GH  COMP9GHOSITE  / 

C..,.  INPUT  MATERIAL  PROPERTIES 

READ(59l000)  LI9L29K9THK, WIDTH 
READ(59 1001)  RHO9UUI29EI9E29GI29GI39G23 
DO  150  J=l93 

DO  150  1=193 

IF(I.EQ.3  .OR.  J.EQo3)  GO  TO  150 
DS(J9l)  =0. 

150  ABD(J9l)  =  ABD(J+3pI)  =  ABDCJpI+O)  =  ABDCJ+OiI+S)  =  0. 

LI  =  MIN0(49MAX0(l9Ll)) 

D(l)  =  LI 

L2  =  MIN0(49MAX0(l9L2)) 

D(2)  =  L2 
D(3)  =  K 
LINT=0 

IF(E1-E2)  12091109120 
no  G12=El/(2o*(l-+UU12)) 

Jl=l  $  J2=3 
GO  TO  200 
120  Jl=4  $  J2=5 

IF(LAYERoEQ.l)  Jl=2  $  J2=3  ^ 

200  WRITECGpEOOO)  LAYER9  WD( JDfWDC J2)*THK9El9E29G129G13»G23»UU129 
RHO9LI9L29K 

CALL  CMPD 
RETURN 

C....  FORMAT  FOR  INPUT-OUTPUT 

1000  FORMAT(3I592F10.0) 

1001  FORMAT(7F10aO) 

2000  F0RMAT(/5X9 129 12H  LAYERCS)  OF92AB92IH  PLATE  WITH  THICKNESS9 

1  FIO.4//IOX9 15HYOUNG?iS  MODULUS9 10X95^El=5^9E10.49l0X95^E2=?^9E10.4/ 

2  10X9 15H3HEAR  M0DULUS9SX9?^G12=?^9  El  0.49  9X9  ?^G13=?^9  El  0.49  9X9 

3  ?^G23=^^9 ElO . 4/10X9 15HP0ISS0N  RATIO9 8X9 ?^UU12=j^9 F5.3/10X9 

4  7HDEN5ITY9  17X9?^RHO=5^9E10. 4/10X9 13HGAUSS  PTS/DIR9 12X9?^L1=?^9 159 

5  5X9  ?^L2=f^9  15/10X9 12HSTRESS  POINT9 14X9  ?^K=?^9 15/) 

END 

C 

SUBROUTINE  CMPD 

COMPUTE  j^ADD?^  MATRIX  AND  ^DS^  MATRIX 
COMMON  /CTDATA/  O9 HEAD(20) 9 NUMNPp NUMEL9 LAYER9 NEQ? IPR 
COMMON  /MTDATA/  RHO9 UU129 El 9  EE? Gi29 G139 G239 THKp WIDTH 
COMMON  /COMPST/  ABDCBp G) 9 DSCEp 2) 9 QBR(39 39 25) 9 QBSCEp 29 25) 9 

-  TH(25)9ZK(25) 

DIMENSION  0(393), QS(292)9TK(25) 

LL=LAYER 

(1M=LL-i-l 

READCSp  1000)  (LpTHCDpTKCDpI^nLL) 

ZKC1)=TTK=0.0 
DO  15  I=l9LL 
TT{<=TTK+TK(I) 

Zi<(T+l)=TK(I)+ZK(I) 

15  CONTINUE 
BO  23  1  =  1 9  MM 
ZK(I)=ZK(I)-TTK/2. 

25  CONTINUE 


PROF  50 
PROF  51 
PROF  52 
PROF  53 
PROF  54 
PROF  55 
PROF  56 

MATL  1 
MATL  2 
MATL  3 
MATL  4 
MATL  5 
MATL  G 
MATL  7 
MATL  8 
MATL  9 
MATL  10 
MATL  11 
MATL  12 
MATL  13 
MATL  14 
MATL  15 
MATL  IG 
MATL  17 
MATL  18 
MATL  19 
MATL  20 
MATL  21 
MATL  22 
MATL  23 
MATL  24 
MATL  25 
MATL  26 
MATL  27 
MATL  28 
MATL  29 
MATL  30 
MATL  31 
MATL  32 
MATL  33 
MATL  34 
MATL  35 
MATL  3G 
MATL  37 
MATL  38 
MATL  39 
MATL  40 
MATL  41 
MATL  42 
MATL  43 

CMPD  1 
CMPD  2 
CMPD  3 
CMPD  4 
CMPD  5 
CMPD  G 
CMPD  7 
CMPD  8 
CMPD  9 
CMPD  10 
CMPD  11 
CMPD  12 
CMPD  13 
CMPD  14 
CMPD  15 
CMPD  16 
CMPD  17 
CMPD  18 


BEL=4.»ATAM(1. )/180. 


DEM  =  1. 

0(1.1)  = 

0(2.2)  =  E2/DEN 

0(1.2)  =  0(2.1)  =  UU12-“Q(2.2) 

0(3.3)  =  G12 

0(1.3)  =  0(2.3)  =  0(3,1)  =  0(3,2)  =  0.0 

OSd.l)  =  G13 

05(2.2)  =  G23 

03(1.2)  =  05(2.1)  =  OoO 

DO  40  1=1. LL 

AMGL=TH(I)5:DEL 

C=COS(ANGL) 

SSa!?d)=Q(l.l)•-C«-=*4■^2.MQ(1.2)+2.■^Q(3.3))*(C*m^^^ 

0BR(1,2. 1)=0BR(2. 1?  I)  =  (Q(1. 1  )+0(2, 2)-4."-fQ(3, 3)  ).f(C^;-W)-''---2 
S  -!-0(1.2)"  (W»«'4  +C"--‘4  )  „  , 

QBR(2. 2. 1  )=Q(  1 . 1  )-"-W“---4+2<.-"-(Q(  1. 2)+2.-»Q(3. 3)  )-~-(G-~-W)^;-^2+Q(2. 2)«C^f---4 
QBR  ( 1 . 3 . 1 )  =OBR  ( 3 . 1 . 1 )  =  ( Q ( 1 . 1 ) -0 ( 1 . 2 )  -2 .  «-Q  ( 3. 3 ) )  C«“-3  + 

$  (Q(1,2)-Q(2.2)-!-2,»Q(3.3))»C«W----3 

082(2.3. 1)=0BR(3. 2. 1)  =  (0(1. 1)--Q(1.2)-2.«D(3.3)  )~-N--- «C-r 
$  (Q(1,2)-Q(2.2)<2.«Q(3.3))-"U-C--3 

QBR(3.3.I)=(Q(l.l)-:Q(2.2)-2o»Q(1.2)-2.«Q(3,3))»(U-»C)«-»2+ 

$  Q(3.3)---(W»-;:-4  +C:;-54  ) 


-  E0«UU12-::-»2/El 
El /DEM 


=  UU12-“Q(2.2) 

=  0(3,1)  =  0(3,2)  =  0.0 


QBS(l.l.I)  =  0S(l.l)»C«"-2  +  QS(2,2)-“W«':f2 
QBS(2.2.I)  =  03(1.  l)-N-"--2  +  03(2. 2)-»C-“»2 
0B3(1.2.I)  =  QBS(2.1.I)  =  (Q5(  1. 1  )-QS(2. 2)  )«-C"W 
40  COHTIMUE 
DO  50  J=1.3 
DO  50  I<=1,3 

ABDU  ^7k’  )=  ABD(J  .K  )+QBR( J.  K,  I  )»(ZK(  I+l  )-ZK(  I ) ) 
ABD(J+3.K  )=  ABD(J  .K+3)=  ABD( J+3,I<)+QBR( J.K. I)« 

$  (ZK(I+l)««-2-ZK(I)«»2)/2. 


+  QS(2.2)-“W»-»2 
+  QS(2.2)-»C-“»2 
(05(1. 1)-QS(2.2))«-C-:^W 


ABD(J+3.K+3)=  ABD(J+3.K+3)+QBR(J.K, I)«(ZK(I+l)«»3-ZK(I)5f«3)/3. 

50  CONTIMUE 
DO  55  1=1. G 
DO  55  J=l,6 

IF(I.GE.3  .OR.  J.GE.3)  GO  TO  55 
IF(ABS(DS(I. J))  .LT.  l.E-OG)  D5(I,J)=0.0 
55  IF(ABS(ABD(I. J))  .LT.  i.E-OS)  ABD(I.J)=0.0 
WRITE (G, 2001)  ((ABD(I, J). J=1.G), I=1.G) 

DO  GO  J=1.2 
DO  GO  l<=1.2 

BO  DS(J°<)  =’dS(J.K)  +  QBS(J,K. I)«(ZK(I+1)-ZK(I)) 

WRITE(G.2002)  ( (DS(I, J) . J=l. 2) . 1=1. 2) 

1000  FORriAT(I5.F5.0,F10.0)  ^ 

2001  FORMAT(//. IX. lOHABD  MATRIX//B(2X. GE13.4/) ) 

2002  FORMAT(/. 1X.9HDS  MATRIX//2(2X.2E13.4/) ) 

RETURN 
END 

C 

3UBR0UTINE  KMLIB 

ASSEMBLE  GLOBLE  ARRAV 
COMMON  G(l) 

DIMENSION  M(l) 

EQUIUALENCE  (G(l).M(l)) 

COMMON  /ISUIDX/  ISW 

COMMON  /CTDATA/  0. HEAD ( 20 ).NUMMP.NUMEL. LAYER, MEQ. IPR 
COMMON  /LODATA/  MDF.NDM.NEN.MST.NKM 
COMMON  /PARATS/  NPAR(14),MEND 
M1=NEND 

N2=Nl+NST-MST-:aPR 
IF(ISW.LE.2)  NE=N2-.'-MKM“-IPR 
IF(ISW.GT.2)  NE=M2+MEafaPR 
CALL  SETMEM(NE) 

CALL  PZERO(G(NEND).NE-NEND)  _  _ _ _ 

CALL  MAS50 1  ( G  ( NPAR  ( 1 ) ) .  G  ( MPAR  ( 2 ) ) .  M  ( NPAR  ( 3 ) ) ,  G  ( NPAR  (  4 ) ) , 

1  M(NPAR(5)),M(NPAR(G)),G(NPAR(7)),G(NPAR(8)),M(MPARO)), 


CMPD  19 
CM PD  20 
CMPD  Si 
CMPB  22 
CMPD  23 
CMPD  24 
CMPD  25 
CMPD  2G 
CMPB  27 
CMPD  28 
CMPD  2j 
CMPD  30 
CMPD  31 
CMPD  32 
CMPD  33 
CMPD  34 
CMPD  33 
CMPD  33 
CMPD  3? 
CI'1P!J  38 
CMPD  ;:;3 
CilPil  40 
CiiPD  41 
Cf1PD  42 
CMPD  43 
CMPD  44 
CMPD  45 
CMPD  4G 
Ci1PD  47 
CMPD  48 
CMPD  43 
CMPD  SO 
CMPD  51 
CMPD  52 
CMPD  53 
CMPD  54 
CMPD  55 
CMPB  S3 
CMPD  57 
CMPD  58 
CMPD  53 
CMPD  GO 
CMPD  El 
CMPD  B2 
CMPD  E3 
CMPD  G4 
CMPD  G5 
CMPD  GG 
CMPD  G7 
CMPD  GS 
CMPD  G3 
CMPD  70 

KMLI  1 
KMLI  2 
KMLI  3 
KMLI  4 
KMLI  5 
KMLI  G 
KMLI  7 
KMLI  B 
KMLI  9 
KMLI  10 
KMLI  11 
KMLI  12 
KMLI  13 
KMLI  14 
KMLI  15 
KMLI  IB 
KMLI  17 
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2  G(NPAR(ll)).G(Nl),G(N2),raF,MDM,NEN,N5T,MKri) 
RETURN 
END 


SUBROUTINE  MASSOKUL^XL^LD^P?  IXp  IDj XpF?  JDIAG^B^Sf  AyNDF»NDM»NEN> 
-  NST,NKM) 

FORM  MASS  MATRIX 

COMMON  /CTDATA/  OfHEAD(20)>NUMNP»NUMEL9 LAVER fNEQ» IPR 

COMMON  /MTDATA/  RH0,UU12,E1»E2.G12.G13,G23,THK, WIDTH 

COMMON  /DMATIX/  D( 10) » DBCB, B), LINT 

COMMON  /ELDATA/  N^NEL^MCT 

COMMON  /ISWIDX/  ISW 

COMMON  /GAUSSP/  SG(1B)»TG(1B),WG(1B) 

DIMENSION  UL(1)»XL(NDM» 1)>LD(NDFj l)»P(l)f IXCNEN? 1)> IDCNDFf 1)» 

1  XCNDM^r  1)9F(1)»  JDIAG(1)»B(1)>S(NST>  l)5»A(l)fSHP(3>  12) 

C....  LOOP  ON  ELEMENTS 
DO  110  N=1,NUMEL 
DO  10  1=1. NST 
DO  10  J=1.N5T 
10  S(I.J)=0, 

C _  SET  UP  LOCAL  ARRAYS 

CALL  PFORMCUL.XL.LD. IX. ID. X. F. B. NDF. NDM. MEN. ISW) 

C..-0  COMPUTE  CONSISTENT  MASS  MATRIX 
L  =  D(l) 

CALL  PGAUSS(L.LINT) 

DO  500  L=1.LINT 

C  . .  COMPUTE  SHAPE  FUNCTIONS 

CALL  SHAPE ( SG ( L ) . TG ( L ) . XL. SHP . XS J. NDM. NEL. IX. .FALSE . ) 

DU  =  WG(L)«-XSJ‘?fRHO-»THK 

C  ..  FOR  EACH  NODE  J  COMPUTE  DB=RHO^»SHAPE»DU 
K1  =  1 

DO  500  J=1.NEL 
Wll  =  SHP(3.J)*DU 
W33  =  Wll'^THl<'^^f2/12. 

C  ..  FOR  EACH  NODE  K  COMPUTE  MASS  MATRIX  (UPPER  TRIANGULAR  PART) 
J1  =  K1 

DO  510  K=J.NEL 

S(J1  .K1  )  =  S(J1  .K1  )  +  SHP(3.K)*W11 

S(Jl+3.Kl+3)  =  S(Jl+3.Kl+3)  +  SHP(3.  K)-»W33 
510  J1  =  J1  +  NDF 
500  K1  =  K1  +  NDF 

C  o.  COMPUTE  MISSING  PARTS  AND  LOWER  PART  BY  SYMMETRY 
NSL  =  NEL^^'NDF 
DO  530  K=1.NSL.NDF 
DO  520  J=K9NSL9NDF 

S(J+29K+2)  =  S(J+1.K+1)  =  S(J  .K  ) 

S(J+49K+4)  =  5(J+39K+3) 

SCK  .J  )  =  S(J  .K  ) 

SCK+S.J'^S)  =  S(J+39K+3) 

S(K+29J+2)  =  S(K+l9J+l)  =  S(J  9K  ) 

520  S(I<+49J+4)  =  S(J+39K+3) 

530  CONTINUE 

IF(ISWoEQo2)  GO  TO  100 
Coc  LUMPED  MAGS  MATRIX 
SUMl  =  OoO 
SUii2  =  OoO 


SUMDl  =0.0 

SUMD2  =0.0 

DO  540  1=1 9 NSL. NDF 

SUMDl  =  SUMDl  +  S(I.I) 

SUMD2  =  SUMD2  +  3(1+3. 1+3) 

DO  540  J=l9NSL.NDF 

SUMl  =  SUMl  +  S(I.J) 

540  SUM2  =  SUM2  +  SCI+O.J+O) 

DO  550  1=1 9 NSL. NDF 


PCI)  = 
P(I+2) 
P(I+3) 
550  P(I+4) 
C. • • «  ADD 


SCI. I) -SUMl /SUMDl 
=  P(I+1)  =  PCI) 

=  S(I+39l+3)^SUM2/SUMD2 
=  P(I+3) 

TO  TOTAL  ARRAY 
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KMLI  IS 
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MASS  1 
MASS  2 
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MASS  4 
MASS  5 
MASS  B 
MASS  7 
MASS  8 
MASS  3 
MASS  10 
MASS  11 
MASS  12 
MASS  13 
MASS  14 
MASS  15 
MASS  IB 
MASS  17 
MASS  18 
MASS  19 
MASS  20 
MASS  21 
MASS  22 
MASS  23 
MASS  24 
MASS  25 
MASS  2B 
MASS  27 
MASS  28 
MASS  29 
MASS  30 
MASS  31 
MASS  32 
MASS  33 
MASS  34 
MASS  35 
MASS  3S 
MASS  37 
MASS  38 
MASS  39 
MASS  40 
MASS  41 
MASS  42 
MASS  43 
MASS  44 
MASS  45 
MASS  4G 
MASS  47 
MASS  48 
MASS  43 
MASS  50 
MASS  51 
MASS  52 
MASS  53 
MASS  54 
MASS  55 
MASS  56 
MASS  57 
MASS  58 
MASS  53 
MASS  GO 
MASS  Bi 
MASS  G2 
MASS  63 
MASS  B4 
MASS  65 
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100  CALL  ADDSTF(A.5»P.JDIAG.LD»tiST»ISEL*NDF»  .FALSE.) 
110  COMTIMUE 

rewind  2  ,  . 

IFdSW.EQ.E)  WRITECE)  (A(I).  I=1.NKM) 
IFdSW.EQ.S)  WRITE(2)  ( Ad ) .  1=1. NEQ) 

RETURN 

END 


SUBROUTINE  RODIPCT 


LOGICAL  FLAG 
COMMON  GCl) 

DIMENSION  MCI) 

COMMON*"/CTDAT  A/  ^  0  ’  HEAD  (20).  NUMNP .  NUMEL .  LA  YER.  NEQ .  I  PR 
COMMON  /LODATA/  NDF.NDM.NEN.NST.NKM 
COMMON  /PARATS/  NPAR( 14) . NEND 
COMMON  /RODATA/  UR. IQ.NDS 

COMMON  /ROELEM/  NER.NEQR.ER  _ 

DATA  FLAG/ . FALSE. /. NER/EO/. ER/30000000 . / 

IF(FLAG)  GO  TO  50 
NEQR=2«(NER+1) 

Ni<MR=7-::-NER+3 

N1=NEND 

N2=N1+NEQ«'IPR 

N3=NE'rNEQ-“IFR 

N4=N3+NEQ«IPR 

H5=N4+N1(MRv:-IPR 

NG=N5+NEQR-IPR 

N7=NS+NEQR 

NB=N7+NEQR»IPR 

N9=N3+NEQR*IPR 

N10=N9+NEQR«IPR 

Nll=N10+NEQR-~-IPR 

NE=N11+NEQR»IPR 

CALL  SETMEM(NE) 

CALL  PZERO (G( NEND ).NE-NEND) 

SO  CALL  UIMPCT(G(NPAR( 1 ) ). G(NPAR(2) ).M(NPAR(3) ).G(NPAR(4) ). 

P  M(NpAR(9)).G(NPARd0)).G(NPAR(ll)).G(Nl).G(N2)» 

1  G(N3).G(N4).G(N5).M(NB).G(N7).G(N8).G(N3),G(N10). 

4  G(Nll)) 

RETURN 

END 

^  SUBROUTINE  WIMPCTCUL.XL.LD.P. IX. ID.X.F. JDIAG.DR.U.B.U.A.RK.RM. 

^  JDR.RU.RUsRA.RB.FR) 

SOLUE  IMPACT  PROBLEM 
LOGICAL  FLAG. TAN 
COMMON  G(l) 

DIMENSION  Md) 

COMMOn''/CTDATA/^o’. Head (20). NUMNP. NUMEL. LAYER, NEQ.  IPR 

COMMON  /TMDATA/  TIME. DT? DDT. FORCE. ALPHA 

COMMON  /LODATA/  NDF.NDM.NEN.NST.NKM 

COMMON  /NITERS/  ITR 

COMMON  /PARATS/  NPAR(14).NEND 

COMMON  /RODATA/  UR, IQ.NDS 

COMMON  /ROELEM/  NER.NEQR.ER 

COMMON  /CONSTS/  AO. A2. A4. AS, AG. A7. AS. AREA 

COMMON  /PROLOD/  PROP 

COMMON  /ISWIDX/  ISW 

KNSION^UL(l),XUl).LD(l),P(l).IX(l).ID(l).Xd),F(l),JDIAG(l). 

1  I]R(1),U(1).B(1).U(1),A(1).RK(1).RN(1)»'JI^I^(1^'^W(1). 

2  RU(1),RA(1),RB(1),FR(1).0(3),QP(3) 

DATA  ITR/5/. FLAG/ . FALSE . /, HIL/l . 4/. INTE/24/ 

if(flag)  go  to  50 

DO  1  1=1,3 
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G7 

riAss 

GS 
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68 
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70 
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7.1 
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72 
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73 
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1 

ROD! 

2 

RODI 

3 

RODI 

4 

RODI 

0 

RODI 

G 

RODI 

7 

ROD! 

8 

RODI 

0 

RODI 

10 

RODI 

RODI 

1 

RODI 

•1  o 

RGDI 

14 

RODI 

15 

RODI 

iB 

RODI 

17 

RODI 

iS 

RODI 

IS 

RODI 

SO 

RODI 

2i 

RODI 

22 

ROD! 

S3 

RODI 

S4 

RGDI 

S5 

RODI 

CIZ- 

RODI 

27 

RODI 

SO 

ROD! 

as 
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30 

RODI 

31 

RODI 

32 

RODI 

33 

RODI 

34 

RODI 

35 

RODI 

3b 

RODI 

37 

WIMP 

i 

WIMP 

2 

WIMP 

3 

WIMP 

4 

WIMP 

5 

HIMP  B 
WIMP  7 
WIMP  8 
WIMP  9 
WIMP  10 
WIMP  11 
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WIMP  13 
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WIMP  17 
WIMP  18 
WIMP  19 
WIMP  20 
WIMP  21 
HIMP  22 
WIMP  23 
WIMP  24 
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QCD^O.O 
QP(I)=0.0 
1  CONTINUE 
IIJ3=1 

TAN=. FALSE. 

REWIND  2 

READ(2)  (B(I),I=1»NEQ) 

FORCE=0.0 

ALPHA=0.0 

PROP=0.0 

NNEQ=NDF*NUMNP 

AO=B./(WIL^^DT)^t*2 

A2=G./(WIL>^DT) 

A4=AO/WIL 

A5=“A2/WIL 

AG=l.-3./WIL 

A7=DT/2. 

A8=DDT/Go 

CALL  FORMROD(RK,RM, JDR) 

DO  10  I=1,NEQR 
10  RUCI)=“UR 
Q(2)=-UR 
FLAG=.TRUE. 

50  ISW=5 

IF(IDSoEQ.NDS)  TAN=.TRUE. 

CALL  FSTREA ( UL, XL» LD, P, IXf IBfXfFf JDIftG» DR> Uf NDF, NDM, NEN> NST» NNEQ) 

DO  20  1=1, NEQ 

A(I)=DR(I)/B(I) 

U(I)=U(I)+DT^‘^A(I) 

U(I)=U(I)+DT^'^U(I) 

20  CONTINUE 
QP(1)=U(IQ) 

QP(2)=U(IQ) 

0P(3)=A(ID) 

DO  30  I=1,NEQR 

RB(I)=RN(I)*(A0*RU(I)+A2«RU(I)+2.*RA(I)) 

30  CONTINUE 

RB I Q=RU ( 1 ) +DT^RU ( 1 ) +DDT/3 . *RA  C 1) 

ROT=0. 000001 
ICOU=0 

DO  100  IT=1,ITR 
RUT=RB  1 0+0(3)  -:-DDT/B . 

AF=-RUT-OP(l) 

CALL  RODLOAD(FIQ,AF) 

DO  110  I=l,NEOR 
FR(I)=RBCI) 

110  CONTINUE 

FR  ( 1)  =FR  ( 1 )  +  ( 1 .  -WIL)  ^FORCE+WILv^F  IQ 

CALL  ACTCOLCRK, FR, JDR, NEQR, . FALSE. , . TRUE. ,0) 

0(3)  =  A4-  ( FR  ( 1 )  -RU  ( 1 ) )  +  A5+^RU  ( 1 )  +  AB^^^R A  ( 1 ) 

RUTT=RB  I Q+0  ( 3 )  --DDT/B . 

ROTR= ADS ( ( RUTT-RUT ) /RUTT ) 

IF(ROTR.LT.ROT)  ICOU=l 
IF(ICOUoGT.O)  GO  TO  200 
100  CONTINUE 
200  DO  210  1=1, NEQR 

FR  ( I )  =A4«-  ( FR  ( I )  “RU  ( I ) )  +A5*RU  ( I )  +A6^^'R A  ( I ) 
RU(I)=RU(I)+DT^RU(I)+A8^KFR(I)+2.^RA(I)) 
RU(I)=RU(I)+A7^KFR(I)+RA(I)) 

RA(I)=FR(I) 

210  CONTINUE 
0(1)=RU(1) 

Q(2)=RU(i) 

0(3)=RA(1) 

FORCE=FIQ 

PROP=FORCE 

ALPHA=-'QC1)-QP(1) 

RODFR=RU(  INTE)’^AREA-ER 

WRITE(8, 8001 )  FORCE, ALPHA, RODFR, (Q( I ) , 1=1, 3) 

8001  F0RMAT(GEi2o4) 
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WIMP  31 
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WIMP  34 
WIMP  35 
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WIMP  49 
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WIMP  54 
WIMP  55 
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WIMP  57 
WIMP  58 
WIMP  59 
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WIMP  61 
WIMP  62 
WIMP  63 
WIMP  64 
WIMP  65 
WIMP  66 
WIMP  67 
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WIMP  69 
WIMP  70 
WIMP  71 
WIMP  72 
WIMP  73 
WIMP  74 
WIMP  75 
WIMP  76 
WIMP  77 
WIMP  78 
WIMP  79 
WIMP  80 
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WIMP  82 
WIMP  83 
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WIMP  85 
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WIMP  90 
WIMP  91 
WIMP  92 
WIMP  93 
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IDS=IDS+1 

IFCIDS.GT.NDS)  IDS=1 
TAN=. FALSE. 

RETURM 

END 

C 

SUBROUTINE  FORMROD(RK» RNrJDR) 

FORM  STIFFNESS  AND  MASS  MATRICES  OF  ROD 
COMMON  /RODATA/  UR,IQ.NDS 
COMMON  /ROELEM/  NER.NEOR.ER 
COMMON  /CONSTS/  A0» A2. A4. A5i. AG* A7. A3» AREA 
DIMENSION  RK( 1) . RM( 1 ) j JDR(2) . D(G) 

DATA  RHOR/.0003225/>RL/1.0/ 

DATA  D/.229 .3Sj .43» .489 .50. .625/ 

EL=RL/NER 

PAI=4.«ATAN(1.) 

JDR(1)=1 
JDR(2)=3 
DO  100  I=l9NER 

IF(I.LT.G)  A=PAI»(D(I)/2.)»«2 

IF(I.GE.G)  A=PAI»(D(B)/2.)»«2 

TT=A»ER/30./EL 

J1=2«(I+1)-1 

J2=J1+1 

J1M1=J1-1 

JlM2=Jl-2 

JDR(Jl)=JDR(JlMl)+3 

JDR(J2)=JDR(Jl)+4 

K1=JDR(J1N2) 

K2=JDR(J1M1)-1 

RKCKl  )=RK(K1  )+TT*3G. 

RK(K2  )=RK(K2  )+TT«3.«-EL 

RK ( K2+ 1) =RK ( K2+ 1) +TT«4 . «EL»»2 
RK ( K2+2 ) =RK ( K2+2 ) -TT»3G . 

RK ( K2+3 ) =RK ( K2+3 ) -TT«3 . «EL 
RK ( K2+4 ) =RK ( K2+4 ) +TT«3G . 

RK ( K2+5 ) =RK ( K2+5 ) +TT*3 . ^EL 
RK  ( K2+G )  =RK  ( K2+G )  -1  T«EL«-«2 
RK  ( K2+7 )  =RK  ( K2+7 )  -TT^^3 .  -EL 
RK  ( K2+8 )  =RK  ( K2+8 )  +TT«4  o  «EL-»»2 
TT=RHOR<-A»EL 
L1=2»I-1 

RMCLl  )=RM(L1  )+TT/2. 

RM(Ll+l)=RM(Ll+l)+TT«EL-»-»2/420. 

RM(Ll+2)=RM(Ll+2)+TT/2. 

RM  ( L 1  +3 )  =RM  ( L 1  +3 )  +TT«-EL««2/420 . 

100  CONTINUE 
AREA=A 

DO  20  I=1.NEQR 
J=JDR(I) 

20  RK(J)=RK(J)+AO»RM(I) 

CALL  ACTCOLCRK.RM, JDR.NEQR. .TRUE. . .FALSE. . 0) 
RETURN 
END 
C 

SUBROUTINE  RODLOADCF. AF) 

COMPUTE  CONTACT  LOADING 
LOGICAL  RELD.UNLD.PIL 

COMMON  /TMDATA/  TIME. DT. DDT. FORCE. ALPHA 
COMMON  /EXDATA/  0(4) 

DATA  UNLD/ . FALSE . /, PIL/ . FALSE . /. RELD/ . FALSE. / 

IF(PIL)  GO  TO  10 

AMAX=AMIN=FMAX=0.0 

PIL=.TRUE. 

10  IF(RELD)  GO  TO  50 
IF(UNLD)  GO  TO  20 
F=Q(1)*AF»»1.5 
IFCF.GE. FORCE)  RETURN 
UNLD=.TRUE. 

AMAX=ALPHA 
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FORM  23 
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FORM  34 
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FORM  37 
FORM  38 
FORM  33 
FORM  40 
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FORM  42 
FORM  43 
FORM  44 
FORM  45 
FORM  4S 
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IF(AMAX.GT.Q(2))  UK=FMAX/((1.-Q(3))*AMAX+Q(2)»Q(3))**Q(4) 
I F  ( AM  AX  o  LE .  Q  ( 2 )  )  UK=FM  AX/ AMAX«‘^*^Q  ( 4 ) 

AMIN=Q(3)*J^(AMAX“Q(2)) 

IF(AMIM.LT.Oo)  AMIM^O.O 
20  IF(AFoLE.AMIN)  GO  TO  30 
F=U!<-(AF-AMIh)-'^^Q(4) 

IFCFoLTc FORCE)  RETURN 
RELD'-^^oTRUE. 

RK=FMAX/ ( AMAX- AM I N ) o  5 
50  IF(AF.LE.AMIN)  GO  TO  30 
F=RK-(AF“ AMIN)-- 1.5 
RETURN 
30  F=0o0 
RETURN 
END 
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